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ABSTRACT 


This dissertation develops basic theory and applications of statistical multirate 
signal processing. Specihc tools and terminology for describing multirate systems in 
the time and frequency domains are presented. An optimal multirate estimator is 
derived in both a direct form and recursive form. The recursive form of the optimal 
estimator allows calculation of the relative change in performance when input signals 
are added or removed from the multirate system. The optimal multirate hltering 
problem also is specialized to the case of optimal multirate linear prediction. An 
efficient method for calculating the multirate linear prediction coefficients and error 
variances is developed through the use of the multichannel Levinson recursion and 
generalized triangular UL factorization. Finally, a multirate sequential classiher is 
derived and applied to the problem of target classihcation. It is shown that classiher 
parameters needed for implementing the multirate sequential classiher are the same as 
those for multirate linear prediction. The methods presented in this dissertation are 
useful for multisensor fusion, particularly when the sensors are operating at diherent 
rates. 
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EXECUTIVE SUMMARY 


Sensor fusion is an important component of target detection and classification. 
The ability to exploit correlated information among multiple sensors with different 
operating characteristics can lead to improved detection and classihcation perfor¬ 
mance. This dissertation addresses multirate sequential classihcation through the 
development of a common multirate framework and the multirate linear prediction 
equations. 

In order to lay the foundation for representing multirate signals, key terms 
applicable to multirate systems are dehned. These are the fundamental rate, system 
period, system phase and decimation factors. These terms are used to explicitly 
describe the periodic nature of the multirate system and allow for the appropriate 
selection of values for the periodic components in the multirate system. In addition, 
fundamental building blocks of the multirate system are presented, along with the 
statistical characterizations of multirate signals as they pass through the multirate 
building blocks. In particular, a linear periodically time-varying hlter is presented 
that permits the input and output signals to be sampled at different rates. 

The hrst multirate optimal hlter considered in this dissertation is the optimal 
linear estimator. The explicit direct form that estimates multiple input signals at 
possibly diherent sampling rates is hrst presented followed by a recursive innovations 
form. This innovations form of the optimal estimator separates the direct form op¬ 
timal hlters into modihed optimal hlters and cross hlters. These cross hlters remove 
any information from one signal that is contained in the other signals, in an ordered 
fashion. In essence, a new set of input signals that are mutually orthogonal are de¬ 
rived and used as the inputs to the modihed optimal hlters. Using the innovations 
form of the optimal hlter allows one to calculate the relative change in performance 
in the optimal estimator when input signals are added or removed from the system. 
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The optimal filtering problem is specialized to the case of optimal multirate 
linear prediction. The multirate Normal equations are derived for a multirate system 
with multiple input and output signals which are observed at different sampling rates. 
For a multirate system with a system period K, there are up to K distinct sets of 
prediction coefficients and error covariance matrices that apply in a periodic fashion. 
An efficient method for calculating the multirate linear prediction coefficients and 
error variances is developed through the use of the multichannel Levinson recursion 
and generalized triangular UL factorization. 

Finally, a multirate sequential classifier is derived starting from the basic the¬ 
ory of sequential hypothesis testing. It is shown that classiher parameters needed for 
implementing the multirate sequential classiher are the same as those for multirate 
linear prediction. A multirate sequential classiher is then implemented and tested us¬ 
ing audio hies of a propellor plane and three A-10 jet aircraft. The experiments tested 
the classiher performance in selecting between the propellor plane and jet aircraft as 
certain system parameters are changed. These parameters are the signal-to-noise 
ratios of the observed signal and the length of the training data. In addition, the 
performance of the multirate classiher is compared to that of the single-channel and 
multichannel classihers using similar data. 
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I. INTRODUCTION 


Within the held of signal processing, a growing interest in the area of multirate 
signal processing has developed over the past three decades. This area focuses on 
the processing of systems that contain multiple signals that can occur at different 
sampling rates. Since this area deals with systems that have components operating at 
different sampling rates, processing techniques and descriptions are needed to account 
for any disparity between these sampling rates. Many advantages of multirate signal 
processing have been found, and techniques have been applied to many areas, such as 
telecommunications, digital audio encoding/decoding, speech and image processing, 
and geophysical signal processing [Ref. 1, 2, 3]. 

Most of the research on multirate signal processing, beyond development and 
characterization of the basic building blocks, has centered on hlter bank theory and 
multirate Kalman hltering [Ref. 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. A 
brief discussion on each of these areas is provided below. In the majority of the 
hlter bank research, there is one input signal (or vector of input signals sampled at 
the same rate). Within the hlter bank structure, the signal is separated into diherent 
subbands sampled at diherent rates. Signal processing techniques are applied, and the 
separate signals are then synthesized into one output signal (or vector of output signals 
sampled at the same rate). This area of research has greatly improved the knowledge 
and understanding of single-input single-output (SISO) and multiple-input multiple- 
output (MIMO) hlter bank structures; however, it has not developed methods for 
processing multiple input signals at diherent sampling rates. 

In the area of statistical multirate signal processing, there are a number of 
papers published. One area that these papers have focused on is the characterization 
of periodic random processes. These papers present methods to describe periodic 
random processes in both the time and frequency domains and discuss such key 
concepts as cyclostationarity of periodic random processes. Another signihcant area 
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of research in statistical multirate signal processing is the development of state-space 
multirate Kalman hlters. These Kalman hlters have been applied to a variety of 
systems, such as small aircraft orbit control, helicopter passive ranging and machinery 
control. In addition, multirate Kalman hlter research has been applied to tracking 
and estimation of mobile robots as well as recovery of lost speech packets in speech 
signal processing. 

This dissertation approaches statistical multirate signal processing from the 
view point of optimal hltering (i. e., Wiener hltering) of multiple input signals at 
different sampling rates. Following the investigation of statistical characteristics for 
multirate signals, an optimal estimator for a desired output sequence is developed. 
This development is then generalized to the linear prediction problem, which leads 
to the general form of the optimal multirate Wiener hlter. Finally, this optimal hlter 
is applied to the sequential target classihcation problem. Much of the foundation 
for the multirate signal processing in this dissertation is derived from [Ref. 16, 17]. 
In addition the optimal estimation and prediction problems in this dissertation are 
multirate extensions of analogous single-channel and multichannel concepts that can 
be found in [Ref. 18, 19]. The multirate sequential classihcation algorithm is derived 
from previous work on the single-channel and multichannel classihcation theory and 
algorithms of [Ref. 20, 21, 22]. 

The following paragraphs review some of the research that has been conducted 
in multirate signal processing. This includes multirate theory, multirate hlter bank 
theory, statistical signal characterization and multirate Kalman hltering. It is by no 
means exhaustive, and indeed much work is still being conducted today. 

A. LITERATURE REVIEW 

The interest in multirate signal processing gained much popularity following 
the 1975 IEEE Arden House Workshop for Acoustics, Speech, and Signal Processing 
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[Ref. 16]. Much of the initial research in multirate signal processing has focused 
on how to describe the multirate system, including the descriptions of some funda¬ 
mental bnilding blocks of the mnltirate system and on how to apply mnltirate signal 
processing to hlter bank theory. In 1983, Crochiere and Rabiner published the hrst 
comprehensive book on multirate signal processing [Ref. 16]. Their work presents 
the basics of sampling rate conversion through decimation and expansion, inclnding 
time and frequency characterization of multirate signals. The application that they 
primarily focus upon is improving efficiency in a mnltirate system with large sampling 
rate changes throngh multistage structures. These multistage structures break the 
sampling rate changes into multiple stages, which can allow relaxation of hlter design 
criteria. In addition to multistage applications, Crochiere and Rabiner also describe 
several forms of nniform hlter banks (e.g., the discrete Fourier transform (DFT) hlter 
bank and the nniform single-sideband (SSB) hlter bank) where the signal is sepa¬ 
rated into mnltiple parts and all parts are decimated by the same value. They extend 
the DFT hlter bank to a generalized non-uniform DFT hlter bank, which allows the 
separated signals to be decimated by diherent values. 

Ten years after the book by Crochiere and Rabiner, the next major publication 
in mnltirate signal processing was by Vaidyanathan in 1993 [Ref. 17]. Vaidyanathan’s 
work addresses the basics of mnltirate signal processing bnt provides greatly expanded 
applications in hlter bank theory. Vaidyanathan presents and thoroughly describes 
many types of multirate hlter banks, including maximally decimated hlter banks, 
parannitary perfect reconstruction hlter banks, linear phase perfectly reconstrnction 
quadrature mirror hlter (QMF) banks and cosine modulated hlter banks. His work 
also presents the concept of periodically time-varying hlters, which plays an impor¬ 
tant part in the processing of mnltirate signals and the problem of qnantization error. 
Lastly, Vaidyanathan introduces the relationship between wavelet transform theory 
and multirate signal processing since the wavelet transform inherently uses nonuni¬ 
form decimation in developing its subbands. Both of these books [Ref. 16, 17] have 
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provided the foundation on which much of the research on multirate signal processing 
since then has been built. 

In addition to these two seminal volumes on multirate signal processing, there 
have been many contributions in journal papers and workshop and conference pub¬ 
lications. Shenoy, Burnside and Parks extend previous work on multirate hlter bank 
theory to optimum minimax hlters [Ref. 1]. By using a generalized Fourier analysis, 
they derived a new error criteria for multirate filter design, which can be used to 
design the optimum minimax multirate hlters for a specihc input signal class. In 
addition, Chen and Vaidyanathan [Ref. 4] extend the concept of polyphase hlters to 
describe rational sampling rate alterations. These polyphase hlters are useful in repre¬ 
senting perfect reconstruction properties of mnltidimensional delay-chain systems and 
periodicity properties of decimated periodic signals. Rules for multirate structures 
were also extended by Evans, Bamberger and McClellan [Ref. 5] to multidimensional 
structnres by hnding greatest common subblattices and computing coset vectors. 

Since the multirate system contains signals at diherent sampling rates, there 
is a natnral periodicity that can be seen in the multirate system. Some basic forms 
of linear periodically time-varying (LPTV) hlters were described by Vaidyanathan; 
however, Saadat Mehr and Chen [Ref. 23] present two methods of describing LPTV 
strnctnres, each of which consists of a periodic switch connected to several linear 
time-invariant (LTI) systems. Saadat Mehr and Chen show that these structnres can 
be used to solve a general approximation problem where an LPTV system with period 
p is approximated by an LPTV system with period p. They extend their work on 
LPTV strnctures to address polyphase and alias-component representations of LPTV 
systems [Ref. 24]. They show that in general a hlter bank can be represented by two 
LPTV systems connected in a cascading manner. They also show that a p-channel 
hlter bank with period m can be represented by an mp-channel LTI hlter bank if p 
and m are relatively coprime integers. 
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Other research on time-varying systems has included work by Phoong and 
Vaidyanathan [Ref. 25] on using a polyphase approach to hlter banks to address 
unusual properties not exhibited by LTI hlter banks. This work is extended in [Ref. 
26] to discuss the factorability of linear time-varying lossless hlter banks. In particu¬ 
lar, they show that all degree-one lossless linear time-varying (LTV) systems can be 
decomposed into a time-dependent unitary matrix and a lossless dyadic-based LTV 
system. Vaidyanathan expanded his work on periodic systems to address periodic sys¬ 
tems with allpass and paraunitary properties [Ref. 27]. By providing a state-space 
representation of periodic LTI systems and introducing the concept of reachability 
and observability in terms of multirate periodic systems, Vaidyanathan shows that 
reachability and observability are not related to the system minimality in a simple 
way, unlike traditional state-space linear systems. In [Ref. 28], Gadre and Patney 
address issues associated with aliasing cancelation and perfect reconstruction within 
a vector context, and they dehne a vector multirate system. These dehnitions lead 
to conditions whereby a vector LPTV system can become a time-invariant system. 
Further, Spurbeck and Scharf [Ref. 29] applied spectral factorization techniques to 
the hlter design for periodically correlated time series. 

Another approach to representing periodic systems was introduced by Misra 
[Ref. 30]. Misra shows that many results from linear time invariant theory can be 
extended to periodic systems. A necessary condition for these results is that an 
equivalent time invariant system must be found for the periodic system. Misra pro¬ 
vides a numerically reliable and simple procedure to hnd the equivalent time invariant 
systems, and he shows that a minimal-order generalized state-space description can 
always be found. 

The research cited up to this point has focused on deterministic signal pro¬ 
cessing. In addition, the primary focus has centered on hlter bank theory and how 
to best separate a signal into multiple parts that are sampled at diherent rates per¬ 
form any necessary signal processing and then resynthesizing the processed signals. 
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In many signal processing techniques, however, applications involve signals that are 
random processes. Consequently, some work has been done to appropriately describe 
the statistical characterizations of multirate signal processing techniques. 

Since multirate systems are inherently periodic in nature, signals within the 
multirate system are usually not wide-sense stationary (WSS). However, under ap¬ 
propriate conditions, stationary input signals can exhibit cyclostationary or wide- 
sense cyclostationary (WSCS) characteristics. In [Ref. 31] and [Ref. 32], Gardner 
presents methods on how to spectrally characterize order cyclostationary signals 
and how to exploit spectral redundancies that might exist in cyclostationary signals. 
Related work in characterization of LTV systems was conducted by Akkarakaran and 
Vaidyanathan [Ref. 33] through the use of bifrequency and bispectrum maps. The 
bifrequency map is a two-dimensional Fourier transform used to characterize the LTV 
system while the bispectrum map is a two-dimensional Fourier transform that charac¬ 
terizes non-stationary random processes. In [Ref. 34] and [Ref. 35], Therrien presents 
methods for dehning correlations functions and power spectra for multirate signals. 
This is extended by Therrien to non-stationary random processes in [Ref. 36] and in 
[Ref. 37]. 

Research into cyclic higher-order statistics of multirate signals has been con¬ 
ducted by Napolitano [Ref. 38]. Napolitano uses cyclic higher-order statistics to 
derive the input-output relations for MIMO linear almost-periodically time-variant 
systems excited by cyclostationary inputs. In addition, he presents a sufficient condi¬ 
tion on the sampling rate to prevent aliasing when reconstructing cyclic higher-order 
statistics of continuous signals from sampled signals. Research into the problem of 
avoiding aliasing in the cyclic higher order spectra of a decimated time series is ex¬ 
tended by Izzo and Napolitano in [Ref. 39] to a generalized form that avoids aliasing 
and imaging effects of a time series decimated by a fractional factor. 

Using cyclostationary spectral analysis, Ohno and Sakai [Ref. 6] developed a 
method of deriving multirate optimal biorthogonal FIR filter banks that minimizes 
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the time-averaged mean-square error (TAMSE) after the high-frequency subband is 
removed. In order to properly describe the hlter bank, cyclostationary spectral anal¬ 
ysis is used since the output of the hlter bank is cyclostationary for a wide-sense 
stationary input. 

Other statistical methods are used by Yurdakul and Dundar [Ref. 40] to 
estimate the quantization error of FIR-based multirate systems. This is practical for 
implementation of multirate hlters since the quantization of hlter coefficients leads to 
quantization errors in the output signal of any discrete-time system. 

One technique of statistical multirate signal processing researched by Sathe 
and Vaidayanathan [Ref. 41] is adaptive hltering. They use the statistical proper¬ 
ties of signals in LTV systems to develop an adaptive hlter structure that is useful 
for identihcation of band-limited channels. A matrix form of this adaptive hlter is 
shown to provide better performance in terms of lower error energy than a traditional 
adaptive hlter but suhers from high computational burden. 

Outside of the signal processing literature, the most prominent technique in 
statistical multirate signal processing available is Kalman hltering. Many applications 
of and techniques for implementing multirate Kalman hlters have been investigated. 
Early work on multirate Kalman hltering by Andrisani and Gan [Ref. 7] uses two 
diherent Kalman hlters in parallel. One of the Kalman hlters processes the high rate 
measurement, at a reduced order. The second Kalman hlter processes the residuals 
of the hrst Kalman hlter in conjunction with the low-rate measurement. Andrisani 
and Gan were able to reduce the computational complexity of their multirate Kalman 
hlter by developing a suboptimal version, although a slight performance penalty was 
incurred. 

Much work in Kalman hltering has been advanced by Bor-Sen Ghen, along 
with You-Lin Ghen and Ghin-Wei Lin [Ref. 2, 3, 8, 9]. Through their work, the 
multirate Kalman hlter has been applied to modeling of autoregressive (AR) and 
moving-average (MA) models through interpolation and estimation of the values of 
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the AR or MA stochastic signals. In addition, these researchers have developed meth¬ 
ods for optimal signal reconstruction in noisy hlter banks and developed a multirate 
Kalman hlter that models system and channel noise. In addition to using a multirate 
Kalman reconstruction hlter, they derived a sample-interpolation algorithm useful for 
the recovery of missing speech packets. 

Further work on Kalman hltering was conducted by De Leon, Kober, Krumvieda 
and Thomas [Ref. 10]. By separating signals into subbands and then implementing 
multiple Kalman hlters, a target can be tracked while reducing the computational rate 
and sensitivity to noisy measurements. Ni et al. [Ref. 11] also developed multirate 
Kalman hlters that are less sensitive to noise. Their model combines the statistical 
model of the input signal with a multichannel representation of the subband signal. 
This form of the multirate Kalman hlter provides the minimum variance reconstruc¬ 
tion of the input signal, provided that the input signal is embedded in the state vector. 
Other applications of the multirate Kalman hlter include orbit control of small air¬ 
craft [Ref. 13], helicopter passive ranging [Ref. 42, 43] and C.N.C. machining control 
[Ref. 12]. In addition, Tornero et al. [Ref. 14] combine the multirate Kalman hl¬ 
ter with a multirate LQG controller to be used in self-location and path-tracking in 
mobile robots applications. 

Another signihcant area of research employing multirate Kalman hltering is 
the multiresolution multirate (MRMR) estimation problem. One of the fundamental 
tools of MRMR techniques is the wavelet transform. The goal of MRMR estimation 
is to improve the overall estimation capability by exploiting features at diherent 
resolutions and sampling rates. Early work in this area was conducted by Basseville 
et al. [Ref. 44] on modeling and estimation of multiresolution statistical processes. 
Chou et al. [Ref. 45, 46] extended this work to recursive estimation and the kalman 
hlter. More recently, Cristi and Tummala [Ref. 15] extended the work even further 
by developing a recursive MRMR Kalman hlter. 

In addition to multirate Kalman hltering, some research on multirate estima- 



tion is available. Haddad, Bernstein and Hnang developed rednced order estimators 
[Ref. 47]. These estimators were derived using a system of equations consisting of one 
modified Riccati equations and two modified Lyapunov equations, and each system 
of equation corresponds to a subinterval of the period associated with the multirate 
system. Another method of estimation was introduced by Hong [Ref. 48, 49]. His 
method focuses on using high-rate measurements to make low-rate estimations. By 
using hlter banks and taking advantage of the lowpass hltering effect of the hlter bank, 
improved approximations of the low rate estimation were obtained. In addition, a 
good approximation of the original measurements was obtained through orthogonal 
transformation by using highpass hltering and lowpass hltering. This thesis had its 
beginning in early work by Therrien and his students and there have been several 
publications along the way. In addition, some of the early work discussed here pro¬ 
vided the groundwork for a companion thesis [Ref. 50] that extended the work along 
various other directions, most particularly to two-dimensional signals with applica¬ 
tions to super-resolution image reconstruction. In [Ref. 51], Cristi, Koupatsiaris and 
Therrien present a method of multirate estimation for a two-signal input along with 
a quantitative analysis of reduction in mean-square error. The multirate estimator 
is extended to m-signals by Kuchler and Therrien in [Ref. 52] (preliminary work for 
this dissertation), and a recursive means of hnding the hlter coefficients and error 
variance as signals are added to the multirate system is presented. In addition, Ther¬ 
rien presents a method of linear prediction for multirate systems in [Ref. 53] along 
with a Levinson-type recursion for hnding the prediction parameters (in conjunction 
with this dissertation). Further work by Therrien and Hawes in [Ref. 54] and [Ref. 
55] present a method of least mean squares (LMS) calculation for multirate systems 
with two input signals. Lastly, some techniques used in the one-dimensional multirate 
system have been extended to the two-dimensional multirate system by Scrofani 
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and Therrien in [Ref. 56] and [Ref. 57]. They use multirate estimation algorithms 
to reconstruct a high-resolution image from a set of low resolution images that are 
subpixel translated (super-resolution). 

B. DISSERTATION OUTLINE 

From the literature review, it is apparent that the development of multirate 
hlter banks has matured greatly. These hlter banks, however, apply to SISO systems 
or MIMO systems where the sampling rates of all the input and output signals are 
the same. This dissertation begins by developing a general framework for multirate 
signal processing. The characterization of random processes as they are applied to 
decimators, expanders and hlters are presented. Some LPTV hlter bank realizations 
are discussed and, in particular, an LPTV hlter bank that allows the input and output 
to be observed at diherent sampling rates is presented. Further simplihcation of signal 
presentation is presented through the use of matrix forms and Kronecker products. 
Whereas the majority of statistical multirate signal processing has focused on the 
Kalman hlter, this dissertation approaches statistical multirate signal processing from 
the point of view of the more general multirate Wiener hlter and its extensions. 

The remainder of the dissertation is organized as follows. In Chapter II, the 
basic framework for describing a multirate system is developed. This includes de¬ 
scribing the appropriate indexing schemes necessary to adequately describe diherent 
signals that are observed at diherent sampling rates and how they can be referenced 
to each other. This chapter also describes the diherent components available to 
multirate signal processing systems. In particular, the expander, decimator, linear 
time-invariant (LTI) hlter and linear periodically time-varying (LPTV) hlter are de¬ 
scribed. Further, the correlation equations for these components are presented. In 
Chapter III, the multirate Wiener-Hopf equations for optimal hltering are developed. 
The multirate optimal hlter is presented both in its direct form and in its innovations 
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representation form. Simulations are conducted to analyze the performance of the 
multirate optimal hlter. In Chapter IV, the multirate optimal hlter is focused on the 
problem of linear prediction. A convenient formulation of this problem is presented 
and a novel efficient algorithm for solving the multirate linear prediction coefficients is 
developed. In Chapter V, the signal processing application of sequential classihcation 
is developed for a multirate system. Simulations are conducted to explore factors 
that affect the overall capabilities of multirate sequential classihcation. In addition, 
performance results for a single channel, a multichannel and a multirate system are 
compared. Finally, in Chapter VI, the research and conclusions of this dissertation 
are summarized and areas of further research are suggested. 
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II. ANALYSIS OF MULTIRATE SIGNALS AND 

SYSTEMS 


A. GENERAL CONCEPTS 

Signal processing techniques have been well developed for processing a single¬ 
input and single-output system, or for processing a multiple-input and multiple- 
output system where all the signals are received at the same sampling rate; the latter 
is referred to as multichannel signal processing. In some signal processing situations, 
however, observed signals are sampled at various different rates, and it may be desired 
to process or exploit these signals together to perform various operations, such as esti¬ 
mation, prediction, detection or classihcation. To provide a basis for addressing some 
of these problems, some new theory needs to be developed that extends the theory for 
single-rate signals and systems. The theory developed for multirate statistical signal 
processing should reduce to the single-channel and single-rate multichannel problems 
as special cases. In this chapter, some basic concepts, notation, and terminology 
for describing stochastic multichannel systems and signals are introduced. In other 
words, this chapter presents the mathematical tools to be used in the remainder of 
the research. 

A representation of a general multirate system is shown in Fig. 2.1. From Fig. 
2.1 it can be seen that each signal has an associated sampling rate, Fj (Hz), where 
the letter i represents a particular signal. These sampling rates are generally differ¬ 
ent, although some may be identical. Since these signals, in general, have different 
sampling rates, their observations do not necessarily occur at the same point in con¬ 
tinuous time and their indices m* may not be aligned. For example, suppose a signal 
X is sampled at 3000 Hz and a signal y is sampled at 2000 Hz. If both signals have 
a sample that occurs at t = 0, then the next observation of the signal x would occur 
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Figure 2.1. Multirate System Concept. 


at 0.33 milliseconds and the next observation of y would occur at 0.5 milliseconds. 
Observations of both x and y would not coincide until 1 millisecond after the hrst 
observation of both signals, this concept is illustrated in Fig. 2.2. A diagram such as 
that provided in Fig. 2.2 will be called a sampling pattern. 


= 3000 Hz 

X • • • «... 

H-^^-1-^^-1-*- time 

0 0.5 ms 1 ms 

y m • • . . . 

F = 2000 Hz 

Figure 2.2. Effect of Different Sampling Rates on Observation Times. 

In order to provide a common framework to which all the signals can be refer¬ 
enced, to account for the effects of sampling signals at different rates, a structure is 
dehned here that will be called the fundamental layer. This structure does not neces¬ 
sarily correspond to any physical portion of the multirate system. It is a conceptual 
tool that is used to allow common referencing of all the signals within the multirate 
system. 
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A key element of the fundamental layer is that it has a sampling rate that is 
a multiple of the sampling rates for all signals in the system. This allows all of the 
signals to be described explicitly in terms of fundamental layer parameters. Return¬ 
ing to the example of two observed signals x and y with sampling rates = 3000 
Hz and Fy = 2000 Hz, respectively, the concept of the fundamental sampling rate 
is illustrated. In order to represent both signals in a common framework (the fun¬ 
damental layer) for analysis, each signal is represented at a higher rate which is an 
integer multiple of the actual sampling rate. For efficiency, the lowest such rate that 
can be applied to all signals is used. In the example just described, this sampling rate 
would be 6000 Hz, and is illustrated by the tick marks on the time axis of Fig. 2.2. 
This sampling rate will be called the fundamental rate, F and its inverse T = 1/F 
will be called the system clock rate. The fundamental rate is the minimum sampling 
rate necessary to describe all signals in the multirate system. It is assumed that each 
sampling rate F) is integer-valued; thus the following dehnition can be made. 


Definition 1. The Fundamental Rate, F, is given by 

F = LCM(Fi,F2,...), (2.1) 

where LCM( ) denotes the least common multiple and Fi, F- 2 , ... represent the sam¬ 
pling rates of all the signals in the system. 

Each actual signal in the system can be associated with (a possibly hctitious) 
signal in the fundamental layer which is sampled at the fundamental rate F. Thus 
each signal in the system can be thought of as a down-sampled or decimated version 
of some signal at the fundamental rate. The decimation factor associated with the 
ith gjgj^al is given by 

A'. = F, (2.2) 

^ i 

i.e., it is the ratio of the fundamental sampling rate to the signal sampling rate for 
the signal of interest. While notation such as x[mx] is used to describe a signal at 
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some sampling rate F^, the notation x[n] will be used to describe the corresponding 
signal in the fundamental layer. Thus these two signals are related by 


x[mx 


x[n 


n=Kxmx 


x[Kxmx 


For a system conceptualized as in Fig. 2.1, a global periodicity exists. The 
system period represents the minimum number of time steps at the fundamental rate 
necessary to establish a repetitive sampling pattern for all signals in the multirate 
system. 


Definition 2. The System Period, K, is the minimum common digital period of all 
elements of the digital system. The system period is given by 


K = LCM{K^,K2,...), 


(2.3) 


where Ki, K 2 , ... represent the decimation factors of all the signals in the system. 

For the signal, the number of samples per period is given by 

Mi = (2,4) 

While K represents the number of time steps at the system rate corresponding to one 
period. Mi represents the number of time steps at the signal rate corresponding 
to one period. 

In order to demonstrate the concepts introduced thus far, the following exam¬ 
ple of a two-channel multirate system is provided. Consider a multirate system that 
has two signals x and y with sampling rates = 3000 Hz and Fy = 2000 Hz, respec¬ 
tively. From (2.1), the fundamental rate for this system is F = LCM(3000, 2000) = 
6000 Hz. Then, from (2.2), the signal decimation factors are calculated to be = 2 
and Ky = 3. Finally, from (2.3), the system period is found to be iF = 6. Figure 
2.3 shows the relationships among these terms. It can be seen from Fig. 2.3 that the 
decimation factors and Ky represent the number of time steps, in the fundamental 
layer, between successive observations (samples) of a particular signal. Notice that 
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Figure 2.3. Illustration of Decimation Factors and System Period. 

the sampling pattern repeats every K = 6 time steps at the fundamental rate. This 
also corresponds to = 3 samples at the rate of x and My = 2 samples at the rate 
of y. Thus given a time n, the relative position within the sampling pattern is the 
same for all times n + iK where i is an integer. This relative position will be called 
the system phase of the system (not to be confused with the phase of its Fourier 
transform). The system phase, denoted by k can be determined for any time n by 
writing 

n = mK + fc, 

where m = \_n/K\. Thus the following dehnition can be made. 

Definition 3. The System Phase, k, is given by 

n=k (mod K), (2.5) 

where n is associated time at the fundamental layer, and k is the residual modulo K. 


B. MULTIRATE SIGNAL PROCESSING 

Some common forms of discrete-time signals can be directly related to an 
associated continuous-time signal. This will be shown using the sampling rates of 
the observed signals of a multirate system and the definition of the fundamental rate. 
Given a continuous time signal x{t), the discrete-time signal for a signal sampled at 
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Figure 2.4. Continuous and Multirate Signals: (a) continuous signal, (b) discrete 
signal at fundamental sampling rate, (c) discrete signal at observed sampling rate. 

the rate F^. would be Xi[mxi] = x{mxiTxi), where T^,. = and i = 1,..., N, where 

Fxi _ 

N is the number of signals in the system. Calculating the fundamental rate as F = 
LCM(Fj,^, Fx 2 , • ■ •, Fxj^), the discrete-time signal associated with the fundamental rate 
would be x[n\ = x{nT), where T = = is the system clock rate. All signals are 
referenced to the initial time of f = 0, thus a:[0] = x(0Tx) = a:(0). For clarity, 
parentheses are used to indicate the time dependence for a continuous signal, and 
brackets are used to indicate the time step for a discrete-time signal. 

The relationships among some various signals are depicted in Fig. 2.4. The 
original continuous signal is denoted by x{t). The signal in the fundamental layer is 
X [n ], which may not exist in reality, and the signal which is sampled at the given rate 
Fx^ is then which is a decimated version of x[n\. 
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Some common multirate signal processing operations are now discussed in 


detail. 


1. Decimation 


The concept of decimation, or downsampling, was introduced earlier in Section 
II.A in conjunction with the dehnition of signals in the fundamental layer. In this 
section, decimation is discussed in a more general sense, as it may relate to various 
signals that may exist in a multirate system. Here decimation is discussed from the 
point of view of obtaining one signal within the system from another. 

The process of decimation eliminates data points in a signal vector in order to 
reduce the sampling rate. The number of data points removed is dependent on the 
decimation factor. Given an integer-valued decimation factor, L, only one of every L 
consecutive samples is retained. 


y[my] = x[Lmy\ , for my = ..., —1, 0,1,.... 


( 2 . 6 ) 


The implied sampling rate for y is thus also reduced by a factor of L, i.e., Fy = F^/L. 
The process of decimation is represented in block diagram form in Fig. 2.5, where L 
represents the decimation factor and my is the index of the decimated signal, such 
that my is the set of integers Z = {..., —1, 0,1,...}. The associated index of the 
input signal is = Lmy. 


-* iL - 


Figure 2.5. Decimation. 


It is well known, [Ref. 17, 58], that the frequency spectrum Y{uiy) of the 
decimated signal is related to the frequency spectrum X{ux) of the original signal by 
the formula 



(2.7) 


The effect of decimation, in the frequency domain, is to stretch the frequency band 
as shown in Fig. 2.6. Notice that if the spectrum of the original signal is nonzero for 
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|a;| > —, aliasing will occur due to undersampling. In most cases this is undesirable, 
L 

although in some cases aliasing can be exploited to obtain useful results. To eliminate 

aliasing, the data sequence may be pre-filtered with an (ideal) low pass filter with a 

vr 

cutoff frequency at |a;| = —. 





Figure 2.6. Frequency Spectrum for Decimation 

Notice that given a signal x[n\ at the fundamental rate, it is possible to generate 
L different decimated signals by applying a time shift to the signal at the fundamental 
rate before decimation. These L signals all have the same rate and are defined by 


[m] = x[Lmy + /], Z = 0,1,..., L — 1. 
The generation of these signals is illustrated in Fig. 2.7, 






iL 






Figure 2.7. Decimation with Time Shift 1. 


where z’' represents a time shift of I units at the input rate. When all of the L possible 
decimation signals are formed, this process is known as complete decimation. Figure 
2.8 illustrates complete decimation for a signal x[n\ at the fundamental rate. 

2. Expansion 

The opposite of decimation is the process of expansion. Expansion, or up- 
sampling, increases the sampling rate by inserting zeros between the data points of a 
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Figure 2.8. Complete decimation illustrated for L = 3: (a) discrete signal at funda¬ 
mental sampling rate, (b) discrete signal at decimated sampling rate, no unit shift 
(Z = 0), discrete signal at decimated sampling rate, one unit shift {1 = 1), (d) discrete 
signal at decimated sampling rate, two unit shifts {1 = 2). 
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signal. The sampling rate is increased by a factor of / if / — 1 zeros are inserted after 
every data point. This process is dehned by 


Vhy] 


{ x[my/I] for TUy div / 
0 otherwise 


( 2 . 8 ) 


where the expression div /” signihes that my is divisible by I, i.e., my/I is an 
integer. The process is represented in block diagram form in Fig. 2.9. The freqnency 


x[m^] 


t / 



Figure 2.9. Expansion. 


spectrum of the expanded signal is related to the original signal by 


y(w„) = 


(2.9) 


The effect of expansion in the frequency domain is to compress the frequency band 
as shown in Fig. 2.10. This compression brings additional copies of the original 






Figure 2.10. Frequency spectrum for Expansion. 


spectrum, known as “images”, into the frequency band of interest (—tt, tt). In order 
to produce a rate change without distortion the expander must be followed by an 

71 

(ideal) low-pass hlter with cutoff frequency at |a;| = y. In the time domain, hitering 
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after upsampling has the effect of “filling in” the zeros with “interpolated” values of 
the signal. 


3. Rate Changes 


Any rate change by a rational factor, J/L, can be achieved by expansion, 


low-pass filtering and then decimation of a signal. (It is assumed that I and L 


are coprime, i.e., these integers have no common factors.) The structure for this 
rate change system is shown in Fig. 2.11. The input x[mx] is at a rate F^. After 
expansion, the rate changes to F^I. The low-pass filter which is designed to operate 
at this new rate {F^I) then serves both to eliminate the images due to expansion and 
provide bandlimiting so the decimated signal will not be aliased. In order to avoid 
aliasing, however, the cutoff frequency of the low pass hlter must be less than or equal 



Fy = F,- I/L. 


Observe that the system rate for this simple system is F = F^I and that 
the decimation factors are given by = I and Ky = L. Thus a rate change is 
accomplished by expanding up to the system rate, hltering and then decimating. 



rate 


rate FI rate F, = FI / L 

X } -^ 


Figure 2.11. Rate Change by a Rational Factor I/L. 
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4. Filtering 

Within the framework of multirate signal processing, there are two broad cat¬ 
egories of hlters that are most useful. These hlters are 1) linear time-invariant hlters 
and 2) linear periodically time-varying hlters, and are discussed in the following sub¬ 
subsections. 

a. Linear Time-Invariant Filters 

Linear time-invariant hlters process a signal at some sampling rate 
and produce an output at the same rate. The input x and output y are related by 
convolution as 


CXD OO 

y[m]= h[k]x[m — k]= x[m\h[n — m] (2.10) 

k=—oo m=—OQ 

where h[m] is called the impulse response and has the same rate as x and y. The 
system is causal if an only if h[m] = 0 for m < 0. In this case the lower limit of the 
variable m in (2.10) can be changed from — cx) to 0. 

The system is stable if and only if 


OO 

\hH 

m=—oo 


< OO. 


The system is said to have hnite impulse response (FIR) if h[m] is a hnite-length 
sequence and said to have inhnite impulse response (HR) otherwise. In the case of 
an FIR hlter (2.10) can be reduced to a hnite sum and the system is always stable. 
b. Linear Periodically Time-Varying Filters 
Linear periodically time-varying (LPTV) hlters are important in multi¬ 
rate applications because of the inherent periodic nature of the multirate system. For 
example, in optimum hltering applications (see Chapter IV), the system periodicity 
requires that the hlter coefficients change in a periodic fashion. This is illustrated in 
Fig. 2.12 where the hlter output y[my] is intended to estimate another signal d[my]. 
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Figure 2.12. LPTV Filter 


The sampling patterns for x, y and d are shown in Fig. 2.12 together 
with a time scale labeled ‘n’ that represents the system rate in the fundamental layer. 
The hltering is assumed to be causal so that at any selected time step only values 
of the input occurring at the same time or earlier on the system time scale are used. 
From Fig. 2.12 it can be seen that at time niy = 6 the input sequence x and the output 
sequence y are aligned, however at time rriy = 7 the input sequence lags the output 
sequence. The difference in system clock time exhibited at these points implies there 
will be a difference in correlation between input and output. Therefore, to estimate 
the output properly at time niy, a new set of hlter coefficients different from those 
at time rriy = 7 are needed, although the hlter coefficients are operating on the same 
input sequence. 

Consider the following specihc example. Given a causal LPTV hlter 
with a hlter order of four, the optimal estimate of d at time ruy = 6 is 

^[6] = /i(°)[0]a;[4] + h(°)[l]a:[3] + h^^\2]x[2] + h^^\3]x[l], (2.11) 

where / = 0,1, 2, 3 are the optimal hlter coefficients for estimating d[6]. At the 

next desired estimate time, rUy = 7, no new input data has been observed. Therefore, 
the estimate ^[7] also is 

^[7] = h(2)[0]a;[4] + h(2)[l]a;[3] + h^^^[2]x[2] + h^^^[3]x[l], (2.12) 
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but / = 0,1, 2, 3 is a different set of coefficients since the filter is estimating 

d[7], instead of d[6]. 

At time ruy = 8, a new observation of the input has already occurred. 
Therefore, the hlter equation at this time instance is 

y[8] = h(^)[0]a;[5] + h^^^[l]x[4] + h^^^[2]x[3] + h(^)[3]a;[2]. (2.13) 

Again the coefficients are different from those occurring in (2.11) and (2.12). Finally, 
at time rriy = 9, the input observation occurs at the same time as the output esti¬ 
mate is calculated. This is the same system phase at the time instance of rUy = 6. 
Therefore, the optimal hlter coefficients at time rUy = 9 are the same as those at time 

rriy = 6. 

Although this example uses specihc time steps to show how an LPTV 
hlter functions, a generalized hlter equation can be found. Let x and y be the repre¬ 
sentations of X and y in the fundamental layer, and let [i] be a set of time-varying 
hlter coefficients. The most generic linear hlter can be written as 

OO 

y[n] = ^ h^^\i]x[n — i]. (2.14) 

^=—00 

Observations of x only occur when its argument is a multiple of however. Thus 
the right-hand side of (2.14) can be rewritten as 

OO 

yW\= X] U^'>[K^ma]x[K^[n/K^\-K^rua], (2.15) 

mx=—co 

where [ J represents the hoor (integer part) of its argument, and the term K^- [n/Kx\ , 
on the right-hand side of the equation, is needed to account for fact that x is not 
observed at every time n. The most recent observation of x, in the fundamental 
layer, is described by Kx\- For example, if is 3, then x is observed only 

at n = 0, 3, 6,.... Thus, at n = 5, for instance, the most recent observation of x is at 
time Kx\n/Kx\ = 3[5/3j = 3 ■ 1 or n = 3. Now recall that a; is a decimated version 
of X that satishes 
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If a corresponding decimated version of the filter sequence is defined as 

m ] = m], 

then the hlter equation (2.15) can be expressed as 

OO 

y[n]= X] h^’'^[m:,]x[[n/K^\ (2.16) 

mx=—oo 

Further, observe that the output y[n\ has values valid only at integer multiples of Ky, 
i.e., n = niyKy. Thus (2.16) becomes 

OO 

y[Kymy]= ^ h^’''>[m:^]x[[Kymy/K^\-ru:^], (2.17) 

mx=—oo 

Finally, by noting that y is a decimated version of y, one can write 

(2.18) 

When the hltering is causal, the lower limit on the sum in (2.18) can be changed from 
rUx = —OO to rUx = 0 since the corresponding impulse response terms are zero. 

Several implementations of LPTV hlters using LTI hlters are possible. 
Two common representations known in the literature are shown in Fig. 2.13 and Fig. 
2.14. Figure 2.13 shows an LPTV hlter realized in a rotary switch form. Given a 
system period K, the rotary switch steps through the sequence 0<k<K — lin 
synchronization with the rate of the output sequence and selecting a different LTI 
hlter to process the input. While this form is typically used in the literature to apply 
when the hlter input and output rates are the same, it easily applies to the more 
general case represented by (2.18). 

Figure 2.14 shows an LPTV hlter realized in hlter bank form. Here the 
output y[n\ is equal to the output of where n = k (mod K). The purpose 

of decimation followed immediately by expansion is to remove any non-zero data not 
associated with the proper phase of the system (for a detailed explanation of this 
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Figure 2.13. Rotary Switch Representation of an LPTV Filter 





Figure 2.14. Filter Bank Representation of an LPTV Filter (From [Ref. 17]) 


filter bank see [Ref. 17]). The hlter bank form of Fig. 2.14 is restricted to situations 
where the input and output are sampled at the same rate. 

A generalized LPTV hlter bank was derived that allows for the input 
x[mx] to be observed at a rate and the output ylmy] to be observed at a rate Fy. 
This more general hlter bank is depicted in Fig. 2.15. This hlter bank operates at 
the fundamental rate F, thus the input signal x[mx] needs to be expanded to the 
fundamental rate by and the output needs to be decimated by Ky. In addition 
the delays are shifted by multiples of the output decimation factor Ky and there are 
My sets of hlters where My is the number of samples per period for signal 

y and rUy = ky (mod My). If and Ky are coprime, then My = K^] but this is 
not true in general. Since x[mx\ is observed at a rate of the clock rate, the hlter 


































coefficients for the filters are of the form 

= ho + hiz~^'" + h2Z~^^^ H-h 

where P is the hlter order. 

The LPTV filter bank of Fig. 2.15 requires a clock that operates at 
the fundamental rate F. If two clocks are available that operate at the input and 
output rates, then the LPTV filter bank can be represented as shown in Fig. 2.16. 
In this form Zy^ = z~^^ and z~^ = z~^^. The decimation and expansion factors are 
replaced with My since the right hand side of the filter is operating at the output rate 
Fy. In addition, the hlters can be described by 

H^^y\zx) = ho + hiz~^ + h2Z~‘^ H-h hp_iz~^^~^\ 


C. STATISTICAL REPRESENTATION OF RANDOM 
SIGNALS 

In order to adequately develop methods and models for processing random 
signals, it is important to be able to describe the statistical characteristics of the 
random signals. Since, in most cases the first and second moments are sufficient 
to describe the necessary statistical characteristics, this section is limited to “second 
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clock rate = F 


Figure 2.15. Filter Bank Representation of an LPTV Filter: Single Clock 
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Figure 2.16. Filter Bank Representation of an LPTV Filter: Dual Clock 


moment” properties. It will be seen that even with this limitation, there are signihcant 
issues to be addressed in the case of multirate signals and systems. 

The hrst order moment, or mean, of a random process is dehned as the expec¬ 
tation of the signal: 

= S{x[m]}. (2.19) 

The mean is, in general, a function of time, although for wide sense stationary (WSS) 
signals it is constant. For cases involving multirate signals and systems the mean 
is (in general) a periodic function of time, i.e., iJix[m] = jjix[m + M] where M is the 
period. 

The second order moments of a random process are dehned via the autocorre¬ 
lation function. Let m and m' represent two different values of the time index. Then 
the traditional autocorrelation function is dehned as 

Rx[rn,m'] = S{x[rn]x*[rn']} (2.20) 

where denotes the complex conjugate. The autocorrelation function is used for 
quantifying second order statistical linear dependence between samples of the signal 
at diherent points in time. 
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The autocovariance function is defined as the autocorrelation function with 
the mean removed: 

Cx[m,m'] = S{{x[m] - Hx[m]){x[m'] - Hx[m']Y}. (2.21) 

Frequently, the autocovariance is more useful than the autocorrelation, although the 
two functions are identical when the mean is zero. The autocorrelation and the 
autocovariance functions are related by, 

Rx[m,m'] = Cx[m,m'] +(2.22) 

By dehning the term lag as the difference between the two time instances, 
I = m — m', the correlation and covariance functions can be represented also in terms 
of the time index, m, and the lag, /, 

Rx[Tn; 1] = S{x[m]x*[m — /]}. (2.23) 

This form is referred to as the “time-dependent” autocorrelation function, or the 
“time-lag” autocorrelation function. An analogous dehnition can be made for the 
autocovariance function. 

Two special cases are important for multirate systems. First, if the random 
process is stationary, at least in the wide sense, then the mean of the random process 
is constant and the correlation and covariance functions are dependent on the lag 
only. Thus the dependence on m in (2.23) can be dropped and the autocorrelation 
function can be written as: 


Rx[l] = S{x[m]x*[m — /]}. 

In this case the traditional autocorrelation function Rx[m,m'], dehned in (2.20) de¬ 
pends only on the difference I = m — m'. 

A second case of importance for multirate systems is the cyclostationary case. 
A more complete discussion of cyclostationarity can be found in [Ref. 31, 32], however. 
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it is simply observed here that for a cyclostationary process the time-lag autocorre¬ 
lation function is a periodic function of m, i.e., 

R. [m; /] = R4m + M; /], 3M (2.24) 

where M is the period. This equation will be taken here as the definition of a cyclo¬ 
stationary process. Note that Rxim] 1] is not necessarily periodic in 1. Although for a 
cyclostationary random process the traditional autocorrelation function R^ [m, m'] is 
periodic in both arguments, with the same period M (See [Ref. 31, 32].) 

For two signals sampled at the same rate the traditional cross-correlation 
function is dehned as 

Rxy[m,m'] = S{x[m]y*[m']}, 

where m and m' are two arbitrary values of the index. The cross-correlation can also 
be represented in the time-lag form as 

Rxy[m] 1] = S{x[m]y*[m - /]}, 

where I = m — m'. Further if the signals are jointly wide-sense stationary, then the 
cross-correlation is a function of its lag term only, 

Rxy[l] = £{x[m\y*[m - /]}. 

In multirate signal processing however, one must be able to hnd the cross-correlation 
of signals at different rates. This requires more careful consideration. The discussion 
here follows the development in [Ref. 35]. 

The cross-correlation function for signals in the fundamental layer can be de¬ 
hned in two forms: 

Rxy\Rxi^y\ ^'[2^[^a:]|/ [^?;]} (2.25) 

and 

Rxy[n-, v] = £{x[n]y*[n - v]}. (2.26) 

The latter is the time-lag form. 
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The correlation function for two signals at different rates can be defined in 
the traditional form. Following the discussion in Section II.B, x[ma] and y[m^ are 
represented in terms of signals in the fundamental layer. This representation makes 
it straightforward to write 

Rxy[mx,my\ = £{x[mx]y*[my\} = £{x[Kxmx\y*[Kymy]} 

Rxy l^R-x^^x: R^y^^y^ • (2.2T) 

One can consider the cross-correlation function Rxy[mx,my] as samples of the 
cross-correlation Rxy in the fundamental layer, as shown in Fig. 2.17. The samples 
are dehned on a rectangular grid or lattice [Ref. 59, 60] of points separated by Kx 
samples of the fundamental layer in the x direction and Ky samples in the y direction 
(see Fig. 2.17). Correlation for values of Ux and Uy not on this lattice are not dehned. 



Figure 2.17. Cross-correlation Lattice. 

Developing the time-lag form of correlation for two signals at different rates is 
somewhat more complex. Consider the following dehnition, which can be written by 
analogy: 

m, = ... — 1 0 1 

Rxy[m]l] = S{x[m]y*[m-l]] = £{x[Kxm]y*[Ky{m-l)]}, 

/ = 0 , 1 ,... 

(2.28) 


33 









With V as the time delay in the fundamental layer, the cross-correlation can be plotted 
in the (n, v) space as shown in Fig. 2.18. The indices m and I map to points on a 
non-rectangular lattice. Note that correlation values do not exist for all possible 
combinations of n and u because samples of y occur only when its argument is an 
integer multiple of Ky. 



■ one period 


V,l 


Figure 2.18. Time-lag Cross-correlation Lattice. 


The lattice is described by two basis vectors vi and V 2 which are defined in 
the fundamental layer as 


Vi = 


K. 

K, - K, 


and 


V2 = 


0 

K,, 


The indices m and I in R^y [m; /] represent a point on the lattice that is reached by 
taking m steps in the vi direction and I steps in the V 2 direction. The time-lag 
correlation function and the traditional cross-correlation function are related as 


RxylR^j Rxy\'^^ ^] • (2.29) 

When the signals x and y are jointly wide sense stationary, then the cross¬ 
correlation as defined in (2.28) is a function of v only (and not n). Therefore, the 
values of correlation depend only on the distance from the n axis in Fig. 2.18 and 
since the sampling pattern is periodic, the values of correlation repeat periodically. 
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Consider the following example. Suppose the cross-correlation function for 
signals in the fundamental layer is given by 

= m<l), (2.30) 

where Ux represents the time instance for the signal x in the fundamental layer and 
Uy represents the time instance of y in the fundamental layer. Note that the symbol 
Rxy is used to distinguish the cross-correlation of the signals in the fundamental layer 
from Rxy, the cross-correlation of the associated decimated signals. 

From (2.27), the cross-correlation of the decimated signals is 

Rxy[mx,my] = S{x[K^mx\y*[Kymy\} 

Rxy\_R-x^^X') R-y^^y\ 

_ a-lK^mx-Kyrtiyl 


Then from (2.29), the cross-correlation can be represented in the (m, 1) space as 


Rxy[m-, 1] = Rxy[m, m-l] 

_ y-\Kxm-Ky{m,-l)\ 

^ y-\(K,,-Ky)m+Kyl\ 


(2.31) 


Note that the cross-correlation of the multirate signals is a subset of the values of 
cross-correlation of the signals in the fundamental layer. 

In order to further illustrate the relation between the {nix, my) and (m,/) 
spaces, the following example is provided. In this example numerical values are used 
to show how valid data points in the {mx,my) space are transformed into valid data 
points in the (m, /) space. Assuming Kx = 2 and Ky = 3 and using the dehnitions 
that m = mx and I = mx — my. Table 2.1 shows the results of the transformations 
and Fig. 2.19 visually depicts how the points map from one space to the other. 
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Table 2.1. Index transformation of the cross-correlation 


{mx,my) 

(0.1) 

(1.0) 

(1.2) 

(2.1) 

(2.2) 


(0;-l) 

(i;i) 

(i;-i) 

(2; 1) 

(2;0) 

Symbol 

0 

□ 

A 

0 

V 


• <13 




A ^ 


^2 -- A 




Figure 2.19. Index Transformation of the Cross-correlation. 


D. INPUT-OUTPUT CROSS-CORRELATION FOR EX¬ 
PANSION, DECIMATION AND FILTERING 

Having defined the mean, autocorrelation and cross-correlation functions, ex¬ 
pressions can be derived for the relationships between these statistics for the input 
and output for some common operations. In particular, the operations of decimation, 
expansion and linear hltering are considered. In the following, it is assumed that 
the autocorrelation function Rx[mx,rn'J of the input as given by (2.20) is known. 
Equivalently, the time-lag form of the autocorrelation Rj;[mx',lx] of the input as de- 
hned by (2.23) is known. Further, if the signal is wide-sense stationary, the time-lag 
autocorrelation function can be reduced to 

Rx[lx\ = S{x[mx]x*[mx - lx]}- (2.32) 
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1. Decimation 

The operation of decimation is depicted in Fig. 2.5 and defined as (2.6). For 
ease of reference, the mathematical description of decimation is also provided below: 

y[my] = x[Lmy\ for m = ..., —1, 0,1,.... 

Given two arbitrary time instances, ruy and niy, the output autocorrelation 
for the decimator is given by 

Ry[my,my] = S{y[my]y*[my]} 

= S{x[Lmy]x*[Lmy]} 

Ry [rUy, rriy] = R^ [Lrriy, LrUy] (2.33) 

Thus decimation results in decimation of the auto-correlation function. 

By dehning a lag term which is the difference between the two time instances, 
I = rriy — m'y and renaming, m = rUy, the autocorrelation function can be represented 
in the time-lag form as 

Ry[m; 1] = S{y[m]y*[m - /]} 

= S{x[Lm]x*[Lm — LI]} 

or 

Ry[m-, 1] = Rx[Lm; LI] (2.34) 

Further, if the signal is wide-sense stationary, then the autocorrelation function is 
independent of m and can be written as 

Ry[l] = Rx[Ll] (2.35) 

The cross-correlation function for decimation is given by 

Rxy[mx,my] = S{x[mx]y*[my]} 

= S{x[mx]x* [Lrriy]} 
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or 


Rxy\j^X':^^y\ Rx\p^x^ (2.36) 

Thus the original cross-correlation function is decimated in one argument only. 
Figure 2.20 visually depicts the effect of decimation on the cross-correlation function. 
The circled locations of the input autocorrelation are the values preserved in the 
cross-correlation of the decimator. 



Figure 2.20. Relation of the Input Autocorrelation to the Cross-correlation of the 
Decimator. 


Using the time-lag dehnition for the cross-correlation of two signals at different 
sampling rates, the cross-correlation of the decimator can be written as 

^xylpT'X) ^yl \R^x ^J/]} 

= S{x[m^]x*[L{m^ - ly)]} 

= £{x[m^]x*[mx - {Lly - {L - l)m3;)]} 


or 

Rxy\p^x^ Rx\p^x^ ^^y (-^ (2.3T) 

The expression for the time-lag cross-correlation function in terms of the auto¬ 
correlation function is rather cumbersome. The geometric representation in the (m, z/) 
is more clear. The effect of decimation on the cross-correlation in the (m; /) space is 
visually depicted in Fig. 2.21. Again, the circled values of the input autocorrelation 
are preserved in the cross-correlation. 
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Figure 2.21. Relation of the Input Autocorrelation to the Cross-correlation of the 
Decimator in the (m; 1) Space. 


Even when the input is wide-sense stationary, the cross-correlation cannot be 
described by its lag term only, since the lag term on the right-hand side of (2.37) 
requires knowledge of m^,. 


2. Expansion 

Again, for ease of reference, the mathematical description of expansion is pro¬ 
vided below: 


y[m] 


x[m/I] 


m div / 


I 0 otherwise 

The autocorrelation of the output signal is dehned as 


Ry[my,m'y] = 8{y[my]y*[m'y\} 

£{x[my/I]x*[m'y/I]} my^m'y div / 

otherwise 


Therefore, 


Ry[my,m'y] = 


Ra,[my/1], m' /I] my, m' div I 


0 


otherwise 


where my and m' represent two arbitrary values of the time index. 
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The time-lag form of the autocorrelation can be written as 

Ry[my-, ly] = S {y[my]y* [uiy - ly]} 

{ S{x[my/I]x*[{my - ly)/I]} TUy, {my - ly) cliv I 

0 otherwise 

Therefore, 

R^[my/I];ly/I] my,ly div / 

Ry[my] ly] — < 

I 0 otherwise 

In the case of the expander, the output is never wide-sense stationary even 
if the input is WSS, because zeros have been added. Therefore, the autocorrelation 
of the expander cannot be represented by a lag only, and a time reference must be 
stated as well. If the input is WSS, the output is cyclostationary however, because 
the autocorrelation function is periodic in its first argument. 

The cross-correlation for expansion is computed as 


Rxy[mx,my] = S{x[mx]y*[my]} 


{ £{x[mx]x*[my/1]} my div I 
0 otherwise 


Therefore, 


Rxy R^X'i ^^y 


{ Rx[mx,my/I] my div I 
0 otherwise 


Thus the original cross-correlation function is expanded in one argument only. 
Figure 2.22 visually depicts the effect of expansion on the cross-correlation function. 
The circled locations of the expander cross-correlation are the autocorrelation values 
preserved from the input. All other values are zero. 
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Figure 2.22. Relation of the Input Autocorrelation to the Cross-correlation of the 
Expander. 


Using the time-lag dehnition for the cross-correlation of two signals at different 
sampling rates, the cross-correlation of the expander can be written as 


Rxy[mx] ly] = 8{x[m^]y*[m^ - ly]} 

j 8{x[mx]x*[{m:, - ly)/I]} 

lo 


{rrix - ly) div / 

otherwise 


{ 8{x[mx]x*[mx - {ly + (/ - l)mx)/I)]} {rrix - ly) div I 

0 otherwise 


Therefore, 


\ R:^[mx;{ly +{I-l)m:^)/I)]} - ly div I 

\ (2.38) 

I 0 otherwise 

The effect of expansion on the cross-correlation in the (m, 1) space is visually 
depicted in Fig. 2.23. Again, the circled values of the cross-correlation are the 
values preserved from the input autocorrelation. The uncircled points represent valid 
locations where the cross-correlation is dehned, but these values are zero. 

Similar to the decimator, when the input of the expander is wide-sense sta¬ 
tionary the output cannot be represented in terms of its lag only. The hrst reason 
is that the lag term on the right-hand side of (2.38) requires knowledge of rrix. In 
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Figure 2.23. Relation of the Input Autocorrelation to the Cross-correlation of the 
Expander in the (m, /) Space. 


addition, the expander inserts zeros, therefore the output is not WSS. It is however, 
wide-sense cyclostationary. 

3. Linear Time-invariant Filters 

For the operation of linear time-invariant hltering the output is given by 

OO 

y[m\= h[r]x[m — r\. 

r=—oo 

Given two arbitrary time instances, m and m', the autocorrelation of the hlter output 

is 


Ry[m,m'] = S{y[m]y*[m']} 


= S ■( i h[Ri\x[m — Ri] I I h[Ro]x[m' — R 

\Ri=—oo / \Rq=—oo 


•0 


OO OO 


GEE h[Ri]h* [Ro]x[m — Ri]x*[m' — Rq] 

,Ri=—oo Ro=—oo 


OO OO 


= E E h[Ri]h*[RQ\Ra\m — Ri,m' — Rq 

Ri=—c>o Ro=—oo 

This can be written as a two-dimensional convolution. 


Ry[m, m'] = Rxim, m] * h[m] * h*[m'] 
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where represents convolution with respect to the associated variable. 

By using the relation Rx[m-, 1] = m — /], where I = m — m', the time-lag 

form of the output autocorrelation becomes 

OO OO 

Ry[m-,l]= ^ ^ h[Ri]h*[Ro]Rx[m-Ri,l - {Ri - Ro)] 

Ri=—oo Ro=—oo 

For a stationary process, the time-lag correlation function is independent of the hrst 
argument. Therefore, the result for stationary processes becomes 

OO OO 

h[Ri]h*[Ro]Rx[l-{Ri-Ro)]. 

Ri=—qo Ro=—oo 

This can be put in the well-known form 


Ry[l] = Rx[l]*h[l]*h*[-l], 

where represents convolution. 

The traditional cross-correlation function between the input and the output 
of the hlter is given by 


or 


Rxy[m,m'] = S{x[m]y*[m']} 


= T < x[m] E h[Ro]x[m — Rq 

I \Ro=-oo 

{ OO 

h*[R^]x[m]x*[m' — Rq] 

Rq=-oo 

OO 

= ^ h*[Ro]Rx[m,m' - Ro] 

Ro=—oo 


Rxy[m, m'] = Rxim, m'] * h*[m'] 

Again using the relation Rxy[m]l] = Rxy[m,m — /], where I = m — m', the 
time-lag form of the output cross-correlation function is: 

OO 

Rxy [m; l]= h* [Ro]Rx[m; I + i?o] 

Ro=—oo 
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or 

Rxy[m-, 1] = Rx[m-, 1] * h*[-l]. 

When the input process is stationary, this reduces to 


Rxy[l] = Rx[l]*h*[-l]. 


4. Linear Periodically Time-varying Filters 

When calculating the correlation functions for a LPTV hlter, care must be 
taken to account for the difference in rates between the input and output. Recall 
that the general equation for an LPTV hlter is given by (2.18). 

The traditional cross-correlation function between the input and the output 
of the LPTV hlter is thus 


Rxy[mx,my] = S{x[mx\y*[my]} 

= S lx[m^] ( ^ h^^'^[R4)]x[[myKy/K^\ - Rq] 


kRo=-oo 


sl h*^^^[R^]x[m,]x*[[myKy/K,\-Ro] 


. Rti=—oc 


= h*^^^[Ro]R,[m,,[myKy/K,\-Ro]. 

Ro=—(X) 

which represents a convolution for with R^ along its hrst argument. 

The autocorrelation for the LPTV hlter is given by 

Ry[my,my] = S{y[my]y*[my]} 

= £ h^’^\R,]x[[myKy/K,\-R,]] ( h^^'\Ro]x[[m'yKy/K,\ - Ro] 


oo oo 


V E E h^^^[Ri]h*^^'\Ro]x[[myKy/K,\ - R,]x*[[m'yKy/K,\ - Rq] 

./?! = —OO Ro = — 00 


oo oo 


E E h^^^[R,]h<^'\Ro]R,[[myKy/K,\ - Ri, [m'yKy/K,\ - Rq], 

Ri=—oo Ro=—oo 
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where my = k (mod K^) and m'y = k' (mod Kx). This equation is the two- 
dimensional convolution along indices my and m'y. 

Because of the floor operation used to transform between the my index and the 
m^ index, the time-lag equivalents for the autocorrelation and the cross-correlation 
functions become rather unwieldly. 

A summary of correlation relations is provided in Tables 2.2 through 2.4. 
A summary of relations in the frequency domain is provided in Appendix A. These 
relations were not used in this dissertation, so are not presented here, but are provided 
in the appendix for reference. 

Table 2.2. Summary of Decimation 
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Table 2.3. Summary of Expansion 


A^x]- 


t / 


'[^v] 


II 

R:,[my/I, m'y/I] my, m'y div I 

0 otherwise 

Ry\my',l^ — 

Rx[my/I-, ly/l] my, ly dlV I 

0 otherwise 

II 

j Rx[mx,my/I] my div I 

1 0 otherwise 

H 

II 

Rxirrix, (ly + (I - l)mx)/I] mx - ly div I 

0 otherwise 


Table 2.4. Summary of Filtering 
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E. MATRIX REPRESENTATION 


In order to develop signal processing algorithms in the time domain, it is fre¬ 
quently easier to work with matrix and vector representations of the systems and 
signals. This section therefore develops tools to represent the basic operations of dec¬ 
imation, expansion and linear hltering. The development extends the ideas described 
in [Ref. 61]. 

1. Decimation 

Consider the decimator depicted in Fig. 2.5. The input and output are related 
by y[m\ = x[Lm\. Define the two vectors of input and output samples 

r 1 ^ 

x= a;[0] a;[l] x[LM — 1] 

y= i/[0] y[l] ■■■ y[M-l] . 

Then the decimation operation can be represented as 

y = Dl,mx, 

where 'Dl^m is an M x ML matrix of I’s and O’s used to extract elements of x as 
observations to be placed in y. The matrix Dl,m will be called a decimation matrix. 
For example, if the decimation factor is L = 3 and the size of the desired output 
vector is M = 4, then the associated decimation matrix is 

100 000 000 000 
000 100 000 000 

Da,4 = . (2.39) 

000 000 100 000 

000 000 000 100 
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This can be represented in a more compact form using Kronecker products. In this 
compact form 

Dl,m = Im ® 

where Im is the M x M identity matrix and is a column vector consisting of a one 
followed by L — 1 zeros, 

dL = 1 0 0 ■ ■ ■ 0 • (2.40) 

Appendix B provides a short discussion of Kronecker products. For a detailed analysis 
of Kronecker products and matrix calculus see [Ref. 62], or for a short tutorial paper 
see [Ref. 63]. 

Using this definition and dropping the subscripts for ease of notation, the mean 
of the decimator is 

lUy = Tjy} = T{Dx} = DTjx} = Dma,, (2.41) 

and the output autocorrelation matrix for the decimator is 

= g{yy*^} = T{Dx(Dx)*^} = DT{xx*^}D*^ = DR,,D*^. (2.42) 

Likewise, the matrix form of the cross-correlation is 

R^j, = ^{xy*^} = T{x(Dx)*^} = T{xx*^}D*^ = R,,D*^. (2.43) 

For two separate signals, Xi and X 2 , that are decimated as 

yi[m\ = a;i[Lim] y 2 [m] = X 2 [L 2 m\, 
the cross-correlation matrix is given by 

^Viy2 = Li,Mi^xiX2^*L2,M2- (2.44) 

where the derivation for (2.44) is similar to that for (2.42) above. 

The covariance matrix is defined as 

Cy = T{(y - mj^)(y - rrij^)*^} = £:{(Dx - Dma,)(Dx - Drria,)*^} 

= DTKx - ma,)(x - 
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or 

Cy = DC,D*^. 

The following relation can also be derived easily, 


Cj = Rj - = D (R, - m,my) D*’’, 


(2.45) 


2. Expansion 

Now consider the expander dehned by the input-output relation 

{ x[m/L] m div L 

■ 

0 otherwise 

The matrix form of expansion can be represented as 

y = 

where the expansion matrix U is an LM x M matrix of I’s and O’s used to insert 
L — 1 zeros between samples of x to form the output vector y. This matrix is called 
an expansion matrix. For example, if the expansion factor is L = 3 and the length of 
the input vector is M = 4, then the associated expansion matrix is 

10 0 0 
0 0 0 0 

0 0 0 0 

0 10 0 

0 0 0 0 

0 0 0 0 

U3,4 = 

0 0 10 

0 0 0 0 

0 0 0 0 

0 0 0 1 

0 0 0 0 

0 0 0 0 
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The expansion matrix can be represented as 

Ul,m = Im ® di, 

where Im is the M x M identity matrix and d; is a colnmn vector dehned in (2.40). 

Using this dehnition and following the procednre used for decimation in the 
previous section, the matrix form of the output autocorrelation for the expander is 

R, = T{yy*^} = UR,U*^. 

Likewise, the matrix form of the cross-correlation is 

R,, = T{xy*^} = R,U*^. 

The matrix results are of little use by themselves because there are large blocks of 
zeros as a result of the expansion operation. When expansion is followed by linear 
hltering however, the zeros £11 in so a more meaningful correlation matrix results. 
The cross-correlation of two interpolated signals, 

{ Xi[m/Li] m div Li 
0 otherwise 

{ X 2 [m/L 2 \ m div L 2 
0 otherwise 

is also of interest. The cross-correlation matrix for these two signals is 

^J/lJ/2 ~ ^{yiy2 } ~ Lx,Mi^XxX2^L2,M2' 

The covariance matrix for expansion follows a similar set of transformations. 
The covariance matrix of the expander is 

Cy = UC,U*^, 
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and satisfies the relation 


= U (R^: - ma;m*^) U*^. 

A direct relationship exists between the decimation matrix and the expansion 
matrix. The expansion matrix is the transpose of the decimation matrix, 

Ul,m = Dl,m- 

The proof of this result follows directly from the dehnition 

Ul,m = Im ® di = (Im ® d|^) = 

where property (3) of the Kronecker products was used (as described in Appendix 
B). This can also be seen in the examples of the decimation and expansion matrices, 
(2.39) and (2.46). 

Some combinations of decimation and expansion are also worth mentioning. 
Expansion followed by decimation (using the same factor) is represented by the prod¬ 
uct matrix Dl,mU l,m = D m- Since this pair of operations results in no change 
of the signal, it must follow that 

= 1m- (2.47) 

On the other hand, decimation followed by expansion selects every term of a 
sequence and replaces all other terms with zeros (the length of the sequence is un¬ 
changed). This operation is represented by the product 

= Im ® d^dl). (2.48) 
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This product matrix is an orthogonal projection matrix. The geometric interpretation 
of this operation is a projection of the vector onto a subspace dehned by every 
component of the original vector. 


3. Filters 

The matrix form of a linear hlter with hnite support 

P+ 

y[m\ = h[r\x[m — r\ (2.49) 

r=P- 


can be represented as 

y = H^x, 

where, the term, H, represents the reversal of the matrix H, which is found by 
inverting the order of the elements along both the columns and the rows (see Appendix 
B.)i 

For the case of a linear time-invariant hlter, the matrix H has the form 


h[P_] 0 

h[l + P_] h[P.] 


0 0 
0 0 


H = 


h[P+] h[P+ - 1] 
0 h[P+] 

0 0 

0 0 


0 0 
0 0 


h[P_] 0 

h[l + P_] h[P_] 


0 0 
0 0 


h[p+] h[p+ - 1] 
0 h[P+] 


^The reversal of a P x Q matrix A with elements aij 


is defined as a matrix A with elements 
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where P_ and P+ are the values of the lower and upper limits, respectively, of the 
hlter in (2.49). Using this dehnition, the matrix form of the output autocorrelation 
of the hlter is 

R, = ^{yy*^} = H^R.Hh 

while, the matrix form of the cross-correlation is 

R,, = ^{xy*^} = R,Hh 

If the lower limit of the hlter matrix P_ equals zero, then the hlter is causal, 
since the hlter relies only on past observations of the input data. On the other hand, 
if —P_ is less than zero, the hlter uses future observations of the input. Therefore, 
the hlter is non-causal and cannot be implemented in real-time. For an LPTV hlter 
there are up to K possible sets of hlter coefficients, The appropriate set of 

hlter coefficients is chosen based upon the system phase, as discussed earlier in this 
chapter. 

F. SUMMARY 

In this chapter, the fundamental building blocks and necessary indexing schemes 
for a multirate system were described. In addition, methods for calculating the cor¬ 
relations and power spectral densities for these building blocks were presented. Using 
these building blocks and relations, it is possible to develop the multirate Wiener- 
Hopf equations for optimal hltering of a multirate system. This is the topic of the 
next chapter. 
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III. LINEAR ESTIMATION 


In this section, the methodology for optimal hltering with mnltirate obser¬ 
vations is developed. The basic principles of linear mean-square estimation apply 
to multirate signal processing; however, due to the periodic nature of the multirate 
system, one optimal hltering equation is not sufficient to fully describe the optimal 
hltering problem. By exploiting the system periodicity of the multirate system, a 
set of optimal hltering equations can be derived, one for each “phase” of the hlter. 
The hlter output (Wiener-Hopf equations) can be a single sequence or a vector of 
sequences based upon the input sequences. 

A. SYSTEM EQUATIONS 

To dehne the problem, assume a multirate observation model with a set of M 
wide-sense stationary observation sequences, Xi[mi] ... These sequences 

represent M observations of the signal s[n] subjected to various linear degradations 
(e.g., additive white Gaussian noise (WGN) or linear distortion), shown in Fig. 3.1. 
These observations sequences are to be used to estimate a single desired sequence 
d[n], as shown in Fig. 3.2. Let Fi,F 2 ,... ,Fm represent the sampling rates of the 
observation sequences and let iFi,iF 2 ,... ,Km be the associated decimation factors. 
It will be assumed that the estimate d[n\ is required at the fundamental rate F, so 
the index ‘n’ is used for this sequence. The desired sequence at time n can be written 
as 

d[n] = d[K ■ I + k] 

where K is the system period, I = \_n/K\ and n = k (mod K). The variable k, is the 
sampling phase, as dehned in Ghapter II, with 0 < k < K — 1. 
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Figure 3.1. M-Channel Multirate Observation Model 


Assume that a causal estimate is desired, i.e., the estimate is calculated using 
only the “present” and “past” observations. The optimal estimate is formed by 
summing the output of linear periodically time-varying (LPTV) hlters, one for each 
observed sequence (see Fig. 3.2). Dehne the vector of hlter coefficients 



- 1 ] 


l<i< M 


0<k< K -1 


where Pi is the order of the hlter and denote the vector of observations from the 
channel by 



Xi[mi] Xi[mi-1] 


Xi[mi - Pi+ 1] 


T 


where rrii = \_n/Ki\. Here m* represents the most recent observation time for Xi if 
the causality constraint is to be maintained. Then the estimate at time ‘n’ is given 
by 

M M 

dk[n] = ^ h[n]; where n = k (mod K). (3.1) 

i=l i=l 

The subscript k in dk[n\ is actually redundant since n = k (mod K)] however, this 
notation is used to emphasize that the statistical properties of the estimate are peri¬ 
odically time varying. 
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Figure 3.2. Multirate Optimal Filter (Direct Form) 


The form of the estimation given by (3.1) and depicted in Fig. 3.2 will be 
called the “direct form” realization since the input sequences are directly weighed and 
summed. In the next section, the methods for hnding the optimal hlter parameters 
in the direct form and the corresponding mean-square error variance are developed. 


B. MULTIRATE WIENER-HOPE EQUATIONS 

Having dehned the estimate, it is now possible to dehne the error as 

ek[n] = d[n] - dk[n] (3.2) 

and to hnd the optimal set of hlter coefficients that minimize T{|£a:HP}- The or¬ 
thogonality principle of linear mean-square estimation [Ref. 18, 58] requires that the 
error be orthogonal to the observation vectors, i.e., 

£{^f\n]el[n]} = Q] i = l,2, ...,M. (3.3) 

Substituting (3.1) and (3.2) into (3.3) and taking the expectation yields 

{ / M \ * 'I M 

j | = [n] = 0 
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where the following terms have been dehned: 


flif’M = (3.4) 

RSfN = £{*!*’(3.5) 


This can be rearranged and written as 

M 

bf i= 1 . 2 ....,Af 

i=i 


or equivalently as 


M 

bf = ER 5 ’'Whf. i=l, 2 ,...,Af. ( 3 . 6 ) 

i=i 

Then using the Hermitian symmetry property (3.6) can be represented 

in matrix form as 


■p (^) 
Kii 

p (fc)* 

■•^12 

p (^)* 




1- 

1—1 

■'^12 

R 

■'^22 

p (^)* 

^2M 


hf 

= 

'-d2 


p(fc)r 

^2M 

p (^) 

■ ^MM_ 




z.ik) 


(3.7) 


The associated error variance according to the orthogonality principle is given by 


al = £{d[n]e\[n]} which, upon substitution of (3.1) and (3.2) becomes 


a 


2 

k 


M 


RM - Y. 

i=l 


~(k)*T 
^ di 



( 3 . 8 ) 


Observe that the mean-square error is periodically time varying, i.e., it is dependent 
on the phase k of the hlter. In order to establish a single hgure of merit for the system. 


an average error is dehned. 





(3.9) 



This geometrical average is chosen, since the error, when expressed in decibels (dB) 
becomes the arithmetic average of the al values in dB. For Gaussian signals, this 
relates to average information (entropy). 


C. CALCULATION OF CORRELATION TERMS 

The correlation terms needed in (3.7) and (3.8) can most easily be generated 
by using the linear algebra concepts described in Chapter II. For any time n and 
corresponding index k, the observation sequence can be expressed as 

(3.10) 

where x* is given by 

~ r U 

^i[n]= Xi[n] Xi[n-1] ••• Xi[n - Pi ■ Ki + 1] 

and is an appropriately-dehned decimation matrix. Note that Xj[n] consists of 
all possible values of Xi[n\ (observed or unobserved) and represents just the 

observed values used by the hlter. 

By virtue of (3.10), the correlation matrix of Xj[n] and Xj[n] is given by 

r|‘> = (3.11) 

where is the correlation matrix of Xj[n] and Xj[n] at the fundamental layer 

Ry = T{xi[n]x*^[n]}. 

The cross-correlation vector between x* and d is 

( 3 , 12 ) 
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where f^j is the cross-correlation of between x* and d at the fundamental layer 

Ydi = £{S^i[n]d*[n]}. 

Note that the dependency on the system phase is completely determined by the 
decimation matrices. 


D. FILTER COEFFICIENTS 


The solution to the multirate Wiener-Hopf equations produces the parameters 
for the direct form realization of the hlter. An alternative realization can be developed 
that leads to some useful insights about the role that each input plays in forming the 
overall estimate dk [n]. This alternative realization is a recursive form of the optimal 
hlter equations, and a complete derivation of these equations is contained in Appendix 
C. It is shown there that the hlter coefficients for the direct form can be written as: 

M 

if - - 5: Ei->- 

j=i+l 




dfc) 


(1 < i < M - 1) 



u{k) _ {~{k) _ p(fc)*rr(fc) 

^‘■M ~ \'-dM '^M '^dM 


where the vector and the matrix BT'' are dehned as 


>{k) 


(3.13) 

(3.14) 



i 

<—1 


1 - 

'^di — 


= 



z.{k) 


R(fc) 


fd{i-l)_ 
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and the matrices and are given by 


= i 


i = 1 



■'^12 

p (k) 


fl{k)*T 

■•^12 

■p (k) 

rl2 

p (k) 

■ ■^2(i-l) 


^{k)*T 

p(fc)*r 
^2(i-l) ■ ■ 

1 



-1 


B]:> i<i<M 


Ef) = 


rSi^ i = 1 

- Gp l<i< M 


The hlter coefficient, in (3.13), can then be re-expressed as 

M 

uik) ^ y.>{k) _ TT{k).{k) 


j=i+l 


(3.15) 


where 


and 


u'(fc) _ ^ (~^{k)*T-r{k) 




With these dehnitions it can be seen from (3.15) that the hlter is comprised of 
a modihed optimal hlter, and cross (or prediction) hlters, ii\j\ which allows 

the optimal hlter to be represented in the form shown in Fig. 3.3. This form is 
referred to as the innovations representation because each branch of the realization 
works on just the new information not present in the other existing channels. It can 
be seen that after applying the cross hlters a modihed set of observations {xj ■ ■ ■ 
is obtained. The modihed observation x[ is the residual after predicting Xi using xj, 
1 < J ^ ~ 1 and represents the new information brought in by channel i. The 

primed observations are mutually orthogonal. The modihed optimal hltering terms. 
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h', i = 1, 2,... ,j, represent the optimal filter for estimating d if Xj[mj], j > i is not 
used. 



(LPTV) 

Figure 3.3. Multirate Innovations Representation 


By separating the filter into cross filter terms and a modified optimal filtering term 
for each input signal, it is possible to determine the improvement in error variance 
that would be derived from adding extra signals, without having to recalculate all of 
the filter coefficients. For each new signal added only the prior filtering terms and the 
modihed optimal filtering terms associated with the new signal are calculated. All 
the previous prior filtering and modified optimal filtering terms remain unaffected. 

By substituting the hlter coefficient dehnition of (3.15) into the error variance 
of (3.8) the error variance can be expressed in a recursive form dependent on the 
number of signals observed. Let al ■ be the phase-periodic error variance when only 
j signals are observed. This form of the error variance will be called the innovations 
form of the error variance equation and is expressed as 
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^k,M ~ '-dM '■'■M 
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-dM 
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M 


-EE E E ■ 
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{k)^ 


-H 


(k)* 


*M-2*M-1 


-H 


(fc)=t 


I’M-! 


M 


'(fc). 

*M 


(3.16) 


The derivation for this form is also contained in Appendix C. Though this form of 
the error variance can become unwieldly for systems with a large number of observed 
signals, this equation is useful to show the improvement in performance brought about 
by each additional channel. 


E. PERFORMANCE STUDY 

A typical scenario for multirate estimation is based on the signal model shown 
in Fig. 3.4. The signal of interest is s[n], therefore d[n] = s[n]. The desired signal is 
subject to independent additive noise in each channel before decimation. In particu¬ 
lar, this study considers a system with two observation sequences Xi[m-[\ and 0 : 2 [m 2 ], 
with decimation factors of iFi = 2 and K 2 = 5. 

Two different types of signals were analyzed for comparison: a second order 
AR process and a periodic signal comprised of two sinusoids. Both signals are dehned 
in the fundamental layer by the equations shown in Table 3.1; where the driving term 
w[n] in the AR process is white Gaussian noise (WGN) with mean zero and unit 
variance. 

The associated autocorrelation functions for each of the signals were calculated and 
are provided in Table 3.2. 
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Figure 3.4. Multirate model with multiple observation sequences 
Table 3.1. Signal Examples 


Signal Type 

Signal Equation 

2 nd 

s[n] = 0.9s[n — 1] — 0.8s[n — 2] -|- w[n\ 

2 Sinusoids 

s[n] = 4cos[0.27rn] -|-2 cos[0.027rn] 


The variances of the two noise sequences rji and 772 were chosen such that the 
Signal to Noise Ratio (SNR) assumed the values, —6 dB, —3 dB, —1.7 dB, 0 dB, 1.7 
dB, 3 dB and 6 dB. The corresponding noise variance can be calculated from 

2 _ R'^(O) 

]^g(SNR/10) 

where Rg is the autocorrelation function for s[n]. A block diagram of the optimal 
estimator along with the sampling pattern for the observations is shown in Fig. 3.5. 
Initial simulations were conducted with hlter orders of Pi = P 2 = 5. 

Figure 3.6 shows the theoretical error variances for the second order AR model, 
computed using (3.8) and (3.9). The hgure compares the error variances for using 
the low-rate signal alone, for using the high-rate signal alone, and for using the low- 
and high-rate signals in combination. In this particular experiment, the high-rate 
signal, Xi, was observed in a 0 dB noise environment and the low-rate signal, X 2 was 
observed with a SNR varying over the range —6 to 6 dB. From this hgure, one can 
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Table 3.2. Signal Autocorrelation Functions 


Signal Type 

Autocorrelation Function 

2nd 

R.(0 = < 

^3.7114 ■ (0.8944)' cos(1.0436/ - 0.0646), 1 > 0 

-3.7114 ■ (1.1180)' cos(1.0436/ - 3.0770), 1 < 0 

2 sinusoids 

Rs{l) = 8cos(0.27r/) -I- 2cos(0.027r/), V/ 


^1 =[^ 1 ] 


Xj = [wj] 


H 


(*) 


d,[n\ = s[n\ 


K,=2 


Channel 1 


-o—o—o-H-o— 


Channel 2 


-o- 


-CH 


-O 


P2 


K^=5 


Figure 3.5. Optimal Estimator Block Diagram 


see that the error variance for the low-rate signal, over the SNR range of —6 dB to 
-|-6 dB, is higher than the error variance of the high rate signal observed at 0 dB. 
Individually, the high-rate signal performs better than the low rate signal. When the 
low-rate and high-rate signals are processed together, however, the error variance is 
lower than for either signal individually. At the lowest SNR simulated (—6 dB) only 
a modest improvement was observed. As the SNR of the low rate signal increases, 
however, the improvement in the error variance approaches 2.5 dB. 

Figure 3.7 shows the theoretical error variances for the sinusoidal model. 
Again, the high-rate signal, Xi, was observed in a 0 dB noise environment and the 
low-rate signal, X 2 was observed with a SNR varying over the range —6 dB to -|-6 
dB. Again in this case, the error variance for the low-rate signal at all SNRs tested 
is higher than the error variance of the high-rate signal at 0 dB. When the low-rate 
and high-rate signals are processed together, however, the error variance is lower than 
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2'^'^ Order AR Process: High Rate SNR = 0 dB 



Figure 3.6. Error Variance vs SNR for the 2^^ Order AR Process 

for either signal individually. At the lowest SNR simulated (—6 dB) only a modest 
improvement was observed. As the SNR of the low-rate signal increases, however, the 
improvement in the error variance approaches 6 dB. 

Additional simulations were conducted with hlter orders Pi and P 2 ranging 
from 5 to 30 for an environment where the SNR for both the low rate and high rate 
signals is 6 dB. Selected results are listed in Table 3.3 for the signal with two sinusoids, 
showing the performance that results when using both the high rate and the low rate 
observations together or using only one set of observations at a time. From Table 3.3, 
one can see that the use of both sets of observations results in signihcant improvement 
(3 dB to 5 dB) over using either xi or X 2 separately, although the error associated 
with hltering using only X 2 is large compared to that resulting from using only Xi. 

The effect on order for the AR process could not be studied. Because of the low 
order of the AR process, varying the hlter orders between 5 and 30 has no appreciable 
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Two Sinusoidal Input: High Rate SNR = 0 dB 



Figure 3.7. Error Variance vs SNR for the Periodic Signal 
effect on the error variance. 

F. SUMMARY 

This chapter addressed the problem of optimal hltering of multiple related 
channels of data observed at different sampling rates. The direct form of the Wiener- 
Hopf equations using linear, periodically time-varying hlters was developed, and the 
error variance of the optimal hlter was observed to be periodically time-varying. In 
addition to the direct form, an innovations form of the multirate Wiener-Hopf equa¬ 
tions was presented. This recursive form provides insight into the relative change 
in performance when additional signals are added or removed from the optimal hl¬ 
ter. Using the geometric average of the periodic error variances as a performance 
measure, it was demonstrated that optimal hltering of multiple channels can provide 
improved performance over optimal hltering using one channel, even if the secondary 
channel when used alone has high error variances. In the simulations conducted, the 


67 




Table 3.3. Error Variances for Periodic Signals with 2 sinusoids (All values in dB) 


Filter Orders 

a;i[n] and 0:2 [n] 

Xi[n}\ only 

X 2 [n] only 

p = 5,g = 5 

-0.43 

3.27 

13.98 

P = 30,g = 5 

-13.54 

-12.26 

13.98 

P = 5,g = 30 

-6.54 

3.27 

9.39 

P = 30, g = 30 

-17.32 

-12.26 

9.39 


largest improvements occured when the signal to be estimated consisted of sinusoids 
in noise. When the signal was derived from a second order AR process, the use of 
multiple channels and/or higher order hlters resulted in only a small improvement in 
performance. It was observed that a single channel hlter of order 5 was sufficient to 
estimate the second order AR signal in noise. These examples used only two obser¬ 
vation sequences, but the derivations and experiments presented here can be applied 
to an arbitrary number of cha nn els. 
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IV. LINEAR PREDICTION 


The previous chapter discussed estimation of data from noisy signals using 
present observations and a hnite number of past observations. Applications that 
require the ability to estimate the present value of the observations themselves using 
only past observations (causal hltering), however, exist. A few of the most common 
applications are focused on speech and image coding, target tracking and target 
classihcation. The process of estimating a signal using only its past observations is 
known as linear prediction. Linear prediction is at the core of all linear estimation 
problems, such as the more general Wiener hlter. Furthermore, specihc insights that 
arise from linear prediction, such as the lattice forms of the hlter, can be carried over 
to the more general problems. 

A. SINGLE-CHANNEL AND MULTICHANNEL RE¬ 
VIEW 

Linear prediction has been well researched for the single-channel and multi¬ 
channel cases. At its most intrinsic form, the goal is to estimate the present value of 
a data sequence from a hnite collection of past data. For the single-channel case, the 
prediction equation can be written 

x[n] = —alx[n — 1] — alx[n — 2] — • • ■ — a*px[n — P] 

p 

= - ^aix[n-i] 
i=l 

where x[n] is the predicted value of x[n] using x[n — 1] through x[n — P] and ai 
through ap are the prediction hlter coefficients (written as conjugate negative values 
for later convenience). The value P is referred to as the prediction hlter order. 
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For the multichannel case, the prediction equation generalizes to 
x[n] = —A^^x[n — 1] — A*^^[n — 2] — ■ ■ ■ — Ap^x[n — P] 

p 

= -^Afx[n-i] (4.1) 

i=l 

where the vector x represents the data observed in each of M channels and Ai through 
Ap are the prediction hlter coefficient matrices. This more general case will now be 
discussed. 

The goal of the linear prediction problem is to hnd an optimal solution for 
the coefficients A* in (4.1), such that the expected value of the squared norm of the 
prediction error e[n] = x[n] — x[n] is minimized. The solution to the linear prediction 
problem results in the well known (multichannel) Normal equations 

R[0] R[l] ■■■ R[P] I E 

R[-l] R|0] ... R[P-1 ] A, ^ 0 

R|-P] R[-F+l] ... R|0] Ap 0 

where the multichannel correlation function is dehned as R[i] = T{x[n]x*^[n — i]} 
and E is the prediction error covariance matrix E = Tje [n]e*'^[n\}. 

If the order of the prediction hlter is allowed to grow at each successive ob¬ 
servation, in order to include all past observations at each step, then the calculated 
prediction errors e[n\ will be orthogonal. If a prediction hlter of sufficiently high 
order is used, then the calculated prediction errors will be approximately orthogonal 
{S{e[i\e*'^[j] = 0} for i ^ j). Thus the prediction error hlter can be thought of as 
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a causal whitening filter. This property is important in autoregressive (AR) modeling, 
moving average (MA) modeling and autoregressive moving average (ARMA) modeling 
(see [Ref. 18, 58, 19].) 


B. MULTIRATE LINEAR PREDICTION THEORY 


The method for extending linear prediction to a multirate system will now be 
considered. Recall from Chapter II, that multirate systems are periodic in nature and 
that one can determine an associated system period K. This system period can be 
used to partition the data into blocks with a common sampling structure as illustrated 
in Fig. 4.1 

Channel 1 

Channel 2 

time 


o- 


j 



-A3 

o- 



—R>— 




Figure 4.1. Multirate System Block Structure 


In this example, Channel 1 is observed at half the fundamental rate, i.e., the 
decimation factor is 2, while for the Channel 2 the decimation factor is 3. The 
system period, as defined in Chapter II, for this example is iC = 6. This has been 
highlighted in Fig. 4.1 by the dashed boxes representing blocks of data. When viewed 
at the fundamental rate, the statistics of the multirate system are cyclo-stationary. 
However, when the multirate system is viewed in blocks, the block statistics are 
stationary. 

In the specific linear prediction problem that is considered here, each new data 
point (in any channel) is predicted as it occurs in time. The prediction is based on 
a finite number Q of previous full blocks of data as well as all of the data within the 
current block that occurs earlier than the point being predicted. Specifically, let m 
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represent a time index for the blocks and i = 1, 2 ,..., P be an index for the data 
occnrring within a block, also ordered in time. For the moment, assnme that data 
points from different channels do not occnr at exactly the same time. The prediction 
eqnation for the {i + 1)'^* data point is then written as 


Xi[m] = —a*iXi-i[m] — ■■ ■ — a*^^_^xi[m\ — a*Jx[m — 1] — ... — a*Qx[m — Q] 
i—l Q 

= N “ £ - (l] ■ (4-3) 

p=l q=l 


The variable Xi[m] represents the observed data point for the system block, 

and x[m] is a vector of all the observed data associated with the system block. 

The elements of x[m] are 

xp[m\ 


x[m] 


xp-i[m] 


(4.4) 


xi[m] 

This method of ordering the block data points is for indexing pnrposes only; it does 
not acconnt for data being observed at the same time. The data is ordered according 
to oldest data first with channel 1 being considered the oldest channel for indexing 
purposes. Figure 4.2 shows these variables for the system depicted in Fig. 4.1. 


Channel 1 


Channel 2 


t[w-l] 


x[ffj] 



x^[m-\] x^[m-\] x^[m\ x^[m\ 


Figure 4.2. Multirate System Block Variables 


The prediction error for Xi [m] is given by 

ei[m\ = Xi[m] — Xi[m] 
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or 




Q 


(4.5) 


Ei[m] = Xi[m] + ^ Q.*pXi_p[m] + ^ a*Jx[m - ij] 

p=l q=l 

By applying the principle of orthogonality [Ref. 18], one can determine the 
prediction error variance for each predicted value and the set of equations that the 
prediction coefficients must satisfy. The prediction error variance satisfies 


i-l 


Q 


= S{xi[m]e*[m]} = R^,[0] + ^ ^ (4.6) 

p=l q=l 

where ra,.x[/] is the vector of cross-correlation terms 


:[/] = 


RxiXp [^] 

RxiXp_l [^] 


RxiXl [^] 

Setting the observations orthogonal to the prediction error produces the equations 


0 = S{xj[m — /]e*[m]} 

p Q 

= RxjXi[~l] + ^ ^ Cli,pRxjXi-p [~^] + '^^^i,q^Xj:>^[Q ~ 

p=l q=l 


(4.7) 


Combining these sets of equations for all observations within a data block, and 
using the shorthand notation 

R; = Rx[/] = £ [x[m]x*^[m — /]] , (4.8) 

produces the matrix form of the multirate Normal equations 


Ro 

Ri 

R-q 


Aq 


E 

R-i 

Ro 

■ Rq-i 


Ai 

= 

0 

R-Q 

Ri-q ■ 

Ro 


_Aq_ 


0 
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where 


0 ■■■0 

1 ■■■0 

.-1,1 0 (4.10) 

ap-i,p-2 ■ ■ ■ 1 

Aj 3-P,j flp—l,i ■ ■ ■ ^ 1) 2, . . . , Q) (^‘ H) 

and the error covariance matrix is 

Up X 

0 aj,_, 

0 0 

The matrix elements represented by an ‘ x ’ do not affect the linear prediction problem, 
so their values do not have to be calculated. 

So far, it has been assumed that data points from different channels do not 
occur at the same time. If points do occur at the same time then they must be jointly 
predicted. For instance, for the system depicted in Figure 4.2, data points occur in 
both channels at the first observation time of the data block. The prediction equation 
for Xi[m] is of the form 

xi[m\ = — aj^x[m — q]. 

<7=1 




1 

cip,i 

Ag = a p^2 
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This equation does not change. The linear prediction equation of (4.3) for X 2 [m], 
however, would change from 


X2[m] = -02,1X1 [m] - y^a2^x[m - q] 

q=l 


X2[m] = - ^ a2jjx[m - q]. 
<7=1 


In other words, Xi[m] and X 2 [m] are predicted simultaneously. The prediction coef¬ 
ficient matrix and prediction error matrix of (4.10) and (4.12) for the 5x5 system 
with no joint observations are 


An — 


1 

0 

0 

0 

«5,1 

1 

0 

0 

as,2 

04,1 

1 

0 

as ,3 

04,2 

a3,i 

1 

as,4 

a4,3 

a3,2 

a2,i 


X 

X 

X 

X 

0 

(rj 

X 

X 

X 

0 

0 

a| X 

X 

0 

0 

0 nl 

X 

0 

0 

0 

0 

a? 


(4.13) 


(4,14) 


These matrices would change to 


An — 


1 

0 

0 

0 

0 

as,i 

1 

0 

0 

0 

as,2 

04,1 

1 

0 

0 

as,3 

04,2 

a3,i 

1 

0 

as,4 

a4,3 

a3,2 

0 

1 


(4.15) 
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and 


(j| X X X X 

0 (J 4 X X X 

E = 0 0 (t| X X • (4.16) 

0 0 0 a2i 

0 0 0 (J 42 erf 

because of the joint observations. Notice that the lower right block in Aq has been 
changed to the 2x2 identity matrix while the lower right block E is changed to a 2 x 2 
error covariance matrix which characterizes the joint prediction of Xi and X 2 - The 
equations for solving the Normal equations must be slightly modihed to account for 
any joint observations. In particular, the covariance elements of the joint observations 
can be calculated from 

Q 

eTij = [0] + ^ ^ (y-i,pRxiXi-j, [0] + 'y ) (4-17) 

p 9=1 

where i and j are jointly observed, the summation over p does not include any joint 
observations and af = af^. In addition, the orthogonality equations become 

0 = S{xj[m — /]e*[m]} 

Q 

= RxjXi[~l] + y ^ Cri,pRxjXi-p [~^] + y ^ ~ !■]■ (4-18) 

P 9=1 

The explicit form of the multirate linear prediction problem can most easily 
be seen by using an example. Figure 4.3 shows a two-channel multirate system with 
decimation factors of Ai = 1 and K 2 = 2 . Channel 1 has two samples per period while 
channel 2 has one sample per period. In addition, observations from both channels 
occur simultaneously at every other system clock interval. Thus joint prediction 
of the channels must occur at those times. The block variables Xi[m] used in this 
example have been labeled in Fig. 4.3. Using a hlter order of P = 1 and assuming 
that the block of observations at m = 0 are available, the linear prediction equations 
for predicting Xi[m] can be determined. For a;i[l] the equations would be 
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Channel 1 

Channel 2 


Xi[2] 

^2 [O] Xj [l] Xj [2] 


m = 0 m = \ 


o— 

Xi[0] X3[0] 

r\ * 

^3- 

X, [l] X 3 [1] 

Pi t 




Figure 4.3. Illustration of Decimation for a Linear Prediction System. 


xi[l] = -ai^ia;3[0] - 01 ^ 32^2 [ 0 ] - ai^ 3 a;i[ 0 ] 

The error would then be e:i[l] = xi[l] — Li[l] or 

ei[l] = a;i[l] + 03^1X3 [0] + 03^22^2 [0] + 01,32^1 [0] 

Applying (4.17) with i = j = 1, the prediction error variance af would be 

= Rxi [0] + (ll,lRxiX3 [1] + Oj1,2RxiX2 [1] + Oi,3-Rxi [l] (4.19) 

and the prediction error covariance for observation X 2 [l] and xi[l] is 

<^21 ~ -Rx2Xi[ 0] + ai^ii?a;2X3 [1] + fll,2-Rx2[l] + Ol,3-Rx2Xi [l] (4.20) 

Applying (4.18), the orthogonality equations for the previous observations are 

-Rx3a:i[~l]+ fll, 1 -^X 3 [0]+ai,2-Rz3X2 [0]+CH,3-Ra;3Xi [0] = 0 (4-21) 

-Rx 2 a:i [~l]+fll,l-Rx 2 a :3 [0]+ Ol, 2 -^X 2 [O]+CH,3-Ra;2a:i [O] = 0 (4.22) 

-Rxi [~l]+fll,l-Rxia :3 [0]+ai,2-RziX2 [0]+ Cq, 3 -^X 1 [O] = 0 (4.23) 

In a similar fashion, the linear prediction equations for X 2 [l] can be written. 

The prediction equation, using (4.3), is 

X2[l] = —a2^;^Xi[0] — 02^2^2 [0] — a2,3^l[0] 
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with an error equation, from (4.5), of 


e 2 [ l ] = X 2 [ l ] + 02,1X3 [0] + 02,2X2(0] + 02,3X1(0] 

The prediction error variance and prediction error covariance equations, derived from 


(4.17), are 

<X2 ~ -^ 2 : 2 ( 0 ] + 02 ,li?a; 2 X 3 (1] + 02 , 2 -Rz 2 (l] + 02,3-Rx2Xi (l] (4.24) 

o"i 2 .^ 0 : 1 x 2 (O] T 02,i7?2;j^3;g (1] T 02,2.Rxix2 (1] T 02 , 3 .^x 1 (1] (^•^^) 

The orthogonality equations for the previous observations, using (4.18), are 

-Rx3X2(~1]+ 02,i7?x3 (0]+O2,2-Rx3X2 (0]+O2,3-Rx3Xi (0] = 0 (4.26) 

-Rx2 (~1]+02,1-Rx2X3 (0]+ 02,2-Rx2 (0]+O2,3-Rx2Xi (0] = 0 (4.27) 

-RxiX2 (~l]+02,l-RxiX3 (0]+O2,2-RxiX2 (0]+ 02,3-Rxi (O] = 0 (4.28) 

Finally, the equations for X 3 (l] can be written as 
X 3 (l] = -03,1X2(1] - as,22^1(1] - “ 3 ,12^3(0] - 03,2X2(0] - 03,3X1(0] 


^3(1] = X3(l] + a3,iX2(l] + as,2X1(1] + 03,1X3(0] + 03,2X2(0] + 03,3X1(0] 

and 

0*3 (0] T a 3 ,l 7 ? 2 ;ga ,2 (0] T a 3 , 27 ? 2 ;g 3 ;j^ (0] T 03,l7?2;g (1] T 03 , 27 ? 2 ;g 3,2 (1] T 03,37?2;ga;j^ (1] ( 4 . 29 ) 

There are no prediction covariance equations since there are no observations occurring 
at the same time as xs)!]. In this case however, the number of orthogonality equa¬ 
tions increases to hve, since the error £ 3 ( 1 ] must be orthogonal to all hve previous 
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observations (see Fig. 4.3) These orthogonality equations are: 


■^ 2 : 2 X 3 [0]T Q;3^1-R2;2 [ 0 ] T Q^3,2-Ra:23:i [ 0 ] Tn3^1-R2;2a;3 [ 1 ] T Cl3^2-Ra:2 [l] Tn3^3-R2;2a;i [l] 0 

(4.30) 

■^xia:3[0]T Q^3,l-Ra:ia:2 [O] T *^3,2-^X1 [O] Tn3^1-R2;^3;3 [l] Tn3^2-Ra:ia:2 [l] T *^3,3-Ra:i [l] 0 

( 4 . 31 ) 

Rx2. [~1]+Cl3^li? 2 : 3^2 [—1]+q;3,2-R XzXi [—1]+ a3,i-Rx3 [0]+a3,2-R 2 : 3^2 [ 0 ] “I“fl3,3-Rx3xi [ 0 ] 0 

(4.32) 

■^ 3 : 23:3 [ C^3, 1 -^X 2 [ 1 ] “l“C^3,2.Rz2a:i [ 1] “l“®3,l-^X2a:3 [O] 4“ ®3,2.^X2 [O] 4~®3,3-^X2a:i [O] 0 

(4.33) 

-^xiX3 [ 1] 4“ Cl3,l-^xiX2[ 1] 4“ Cl3,2.^X1 [ 1] 4“ fl3^l/?2;iX3 [O] 4“ fl3,2.Ra;iX2 [0] 4“ ^3,3-^X1 [O] 0 

(4.34) 

By combining (4.19) through (4.34), the matrix form of the multichannel Nor¬ 
mal equations (4.9) can be produced and written as 



where 

Rq Ri 

R = . (4.35) 

R-i Rq 
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or, more explicitly, 


Rx3 [0] 

Rx3X2 [O] 

-^2:3X1 [0] 

: Rx3[1] 

Rx3X2 [1] 

Rx3Xi [1] 

Rx2X3 [O] 

Rx2 [0] 

-Rx2Xi [0] 

■ Rx2X3 [1] 

Rx2[1] 

Rx 2 Xi [ 1 ] 

RxiX^ [O] 

-^0:10:2 [0] 


RxiX3 [1] 

RxiX2 [ 1 ] 

RxA^] 


Rx3X2 [~1] 

Rx3Xl [~1] 

: i? 3,3 [0] 

Rx3X2 [0] 

-Rx 3 Xi [0] 

Rx2X3 [ 1] 


Rx2Xi[ 1] 

Rx2X3 [0] 

Rx2[0] 

-Rx2X1 [0] 

_RxiX3 [~1] 

RxiX2 [ 1] 


-^XiX 3 [0] 

-RxiX 2 [0] 

-RzjO] _ 



1 

0 

0 


0^3,1 

1 

0 

Ao 


«3,2 

0 

1 

Ai 


^3,1 

^2,1 

^1,1 


^3,2 

02,2 

^1,2 


0 - 3,3 

02,3 

Ol ,3 


and the error matrix is 




X 

X 


0 



E 


0 

^?2 


0 


0 

0 

0 


0 

0 

0 


0 

0 

0 


(4.36) 


(4.37) 
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C. AN EFFICIENT SOLUTION TO THE MULTIRATE 
NORMAL EQUATIONS 


In the previous section, the Normal equations for the multirate linear predic¬ 
tion problem were developed element by element. Although these equations can be 
solved directly, it is desirable to develop a method for calculating the linear predic¬ 
tion coefficients and error variances more efficiently. Efficient recursive methods for 
solving these equations have been developed for both the single-channel and mul¬ 
tichannel cases. In the single-channel case the method is called the Levinson (or 
Levinson-Durbin) algorithm while in the multichannel case it is known commonly as 
the multichannel Levinson algorithm, or in some literature it is called the Levinson- 
Wiggins-Robinson (LWR) algorithm [Ref. 64]. In the following, a similar procedure 
is developed for the multirate case. 

The equations to be solved are 


Ro 

Ri 

Rq 


Aq 


E 

R-i 

Ro 

■ Rq-1 


Ai 

= 

0 

R-q 

Ri-q ■ 

Ro 


_Aq_ 


0 


These are the multirate Normal equations of (4.9), repeated here for convenience. 

To begin the solution to (4.38), consider the problem of predicting an entire 
block of data at once. This is identical to the multichannel problem. The Normal 
equations for this problem are of the form (4.2) and can be written as 


- 1 

O 

o 

_1 


I 


1 - 

H 

1_ 

R-1 Ro • ■ ■ Rq-1 


a; 

= 

o •• 

R-Q Ri-Q ■ ■ ■ Ro 




- 1 

• o 


(4.39) 
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Observe that the correlation matrix in (4.39) is identical to that in (4.38) although 
the variables A) and E' are not. If (4.39) is right multiplied with Aq however, and 
(4.38) compared to (4.39), it can be seen that the following equalities must hold: 

A, = A;Ao, ^ = l,...,g (4.40) 

E = E'Aq. (4-41) 

Therefore, if solutions to (4.39) and the matrix Aq can be found, then the desired 
multichannel parameters A* (for i = 1,2,... ,Q) and E can be found from (4.40) and 
( 4 . 41 ). 

A computationally efficient solution to (4.39) can be found using the LWR 
algorithm (see Appendix D). This solution produces the matrices A), where i = 
1,2,... ,Q and E'. To hud Aq from (4.41) the forms of E and Aq must be considered. 
In particular, E is upper triangular, e.g., (4.14), when there are no simultaneous 
observations, and block upper triangular, e.g., (4.16), when some observations occur 
simultaneously. Also, Aq is lower triangular with unit diagonal, e.g., (4.13), when 
there are no simultaneous observations, and block lower triangular in general, e.g., 
(4.15). Consider the specihc case shown in Fig. 4.2. By writing (4.41) as E' = EAg 
the matrices have the form 


/2 

^5 

/2 

*^54 

/2 

<^53 

/2 

<^52 

/2 



X 

X 

X 

X 


1 

0 

0 

0 

0 

/2 

*^45 

a 

^4 

/2 

<^43 

a 

<^42 

a 

*^41 


0 

aj 

X 

X 

X 


0^5,1 

1 

0 

0 

0 

/2 

<^35 

/2 

*^34 

/2 

^3 

/2 

*^32 

a 

*^31 

= 

0 

0 


X 

X 


tt5,2 

«4,1 

1 

0 

0 

/2 

<^25 

a 

<^24 

/2 

<^23 

/2 

t^2 

a 

<^21 


0 

0 

0 




“5,3 

^4,2 

^3,1 

1 

0 

/2 

.^15 

/2 

0 ' i 4 

/2 

<^13 

a 

<^12 

/2 

. 


0 

0 

0 

f^?2 



tt5,4 

^4,3 

0 ^ 3,2 

0 

1 


(4.42) 

It can be seen that (4.42) is similar in form to M = UL, where U and L are upper 
and lower triangular matrices, respectively. Thus if E' can be factored into E and 
Ag then the solution to the multirate Normal equations has been found. 

Most linear algebra books (e.g., [Ref. 65]) contain LU factorization algorithms 
for solving the decomposition of a matrix into a lower-upper triangular scheme (i. e.. 
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M = LU). For the purpose of this dissertation, however, it is desired to factor E' 
into an upper-lower scheme. In addition, this factorization must be general enough to 
allow for block elimination in the algorithm (this accounts for joint observations in the 
multirate system). Thus a generalized UL factorization algorithm has been written 
and provided in Appendix E. Using this generalized UL factorization algorithm to 
hnd Aq, the multirate Normal equations can be solved. A signihcant advantage of 
decomposing E' into E and Ag ^ is that the only matrix that requires inversion is Ag 
Since, Ag ^ is lower triangular, Ag can be found efficiently using forward elimination 
which does not explicitly perform a full matrix inversion [Ref. 66]. 

The steps needed to solve for the multirate Normal equations are summarized 
in Table 4.1. 


Table 4.1. Steps in Solving the Multirate Normal Equations 

1. Solve (4.39) using the multichannel Levinson recursion to hnd E', A'^, ..., Ag 

2. Factor E' using the generalized UL factorization algorithm 
to obtain E and Ag ^ 

3. Invert Ag ^ to obtain Ag 

4. Compute A* = A)Ag for i = 1,..., Q 
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D. SUMMARY 


In this chapter, the multirate linear prediction equations were developed. Also, 
an efficient method for calculating the linear prediction coefficients using the multi¬ 
channel Levinson recursion and LU factorization was propsed. These techniques for 
calculating the linear prediction coefficients are used in the next chapter on multirate 
sequential classihcation. 
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V. CLASSIFICATION 


The problem of distinguishing between two or more different types of signals 
(arising from different sources) is known as classification. The classification is said to 
be “sequential” when the signal to be classified is observed one discrete sample at a 
time. After a new sample is observed, an attempt is made to classify the signal using 
all of the samples available up to that time. If a classification cannot be made, then 
another sample is taken and the procedure continues. The goal is to classify the signal 
with a prescribed probability of error using as few samples as possible. The general 
theory of sequential classification was originally and extensively developed by Wald 
in 1947 [Ref. 20]. It has since been applied specifically to statistical signal processing 
and described for the single-channel and multichannel cases (e.g., [Ref. 22] and [Ref. 
21].) The focus of this chapter is to demonstrate the feasibility of developing an 
algorithm that allows for multirate observations to be used in a sequential classifier. 

A. SEQUENTIAL CLASSIFICATION 

Before developing the sequential classifier for the multirate system, it is nec¬ 
essary to provide a brief description of the sequential classification process. Conven¬ 
tional methods of classification usually involve a fixed block of data. This method 
is known as classifying with a fixed sample size. Statistically optimal methods then 
involve developing the ratio of likelihood functions for the two classes and comparing 
that likelihood ratio to a fixed threshold. This threshold is chosen to optimize some 
criterion, such as total probability of error (Bayes test) or maximizing the probability 
of correctly choosing one class while holding the probability of error on the other 
class fixed (Neyman-Pearson test) [Ref. 22, 67, 68]. Wald developed a method, called 
the sequential probability ratio test (SPRT), which allowed one to successively add 
samples to data set being used to perform the classification. If classification is not 
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possible at that time, another sample is added. This process is repeated until a clas- 
sihcation can be made with some hxed desired probability of error. For an extensive 
discussion on sequential analysis, see [Ref. 20]. 

To begin the discussion of sequential classihcation, assume that the signal of 
interest, x, belongs to one of only two classes, class 0 or class 1. Let Hq be the 
hypothesis that the signal belongs to class 0 and Hi the hypothesis that the signal 
belongs to class 1. Let the vector of observations be denoted by 


x„ = 


^0 




(5.1) 


where Xj is the set of observations of x available at time j, and let /(x„|iLj) be the 
probability density function for x„ given that the observations are from class i. In 
the Wald SPRT, the likelihood ratio is dehned as 


^ /(XnlLTl) 
/(xnli^o) 


(5.2) 


and two positive constants or thresholds A and B (with A > B) are chosen to classify 
the signal x„ as either class 0 or class 1 with some desired probabilities of error. The 


test performed is as follows: if 

f{^n\Hi) ^ ^ 

/(x„|iLo) “ 

then x„ is classihed as class 1; while if 


(5.3) 


f(XnlHo) ^ 


(5.4) 


then Xn is classihed as class 0. If the ratio of the conditional probabilities falls between 


A and B, that is, if 


B < 


/(Xn|Ffl) 

/(x„|iLo) 


< A, 


(5.5) 



then another sample is taken and the process is repeated with 


^n+l ~ 

■hn+1 

Proper selection of the values for A and B is discussed below. 

In many cases, it is more convenient to work with the logarithm of the like¬ 
lihood ratio rather than the likelihood ratio itself. Since the logarithm is a strictly 
monotonic function, this does not change the essential nature of the SPRT. In the 
particular case that follows, it is most convenient to deal with the negative logarithm 
and define 

h„(x„) = -2 ln(/„(x„)) (5.6) 

with 



Ta = -2 ln(R) 

(5.7) 

Tb = -2 ln(R) 

(5.8) 


The inequality then reverses and the SPRT decision becomes 

< ta choose Hi 

h„(x„) 

> tb choose Hq 


where 


hnip^n) — 2 In 


/x„|Hi (XnliPl) 
f:x.„\Ho{^n\HQ) 


Otherwise take another observation. 


(5.9) 

(5.10) 


The choices for the thresholds A and B are related to the probabilities of error 
for class 0 and class 1. Using the method described in [Ref. 22, 21], the relations 
between thresholds and probabilities of error can be derived as follows. Let Cq be the 
probability of a decision error given the data is from class 0 and ei be the probability 
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of error given the data is from class 1. Assume that the upper threshold is crossed 
after n observations; that is, at time n the ratio ln{^n) dehned by (5.2) is greater 
than but approximately equal to A: 


Li^n) > A 

Substituting (5.2) into (5.11) and rearranging yields 

> A/(x„|iJo)- 


(5.11) 


(5.12) 


The condition (5.12) dehnes a region of the x„ space (Ri) where the data is classihed 
as class 1 {Hi). Integrating both sides over this region yields 



/(x„|Tri)dx„ > A 



/(x„|h/'o)dx„- 


(5.13) 


Now, since Ri is taken to be the region of measurement space where x„ is classihed 
as class 1, the left hand side of (5.13) is the probability of correctly classifying class 
1, namely 



/(x„|i^i)(ix„ 


1 — €i. 


( 5 . 14 ) 


The integral on the right hand side of (5.13) is then the probability of incorrectly 
classifying class 0, 



/(x„|Tro)dx„ 


eo 


( 5 . 15 ) 


Substituting (5.14) and (5.15) into (5.13) yields 


1 ~ ^ Aeo 


(5.16) 


or 

-1 

A = (5.17) 

eo 

where the inequality has been replaced by an equal sign since the threshold A is most 
restrictive at that value. Similarly, assuming 


l„(x„) < B 


( 5 . 18 ) 



and integrating over the region where x„ is classihed as class 0 yields 

B = (5.19) 

i — eo 

Eqnations (5.17) and (5.19) provide the necessary relations to allow one to choose 
appropriate values of thresholds A and B based on desired probabilities of error. 

B. ALGORITHM DEVELOPMENT 

When the observations in the SPRT are jointly Gaussian for both classes, 
then a convenient recursive classification algorithm can be developed that involves 
linear prediction. This recursive form has been developed in [Ref. 22, 21] for the 
case of a single observation sequence, which is scalar-valued or vector-valued (the 
multichannel case). Here, the method is extended to the case of multiple channels 
sampled at different rates. Thus, at any given epoch on the system time scale, there 
may be one or more simultaneous observations available from the various channels 
or possibly no observations. The algorithm to be described takes this more general 
situation into account. 

The development of the recursive form of the quadratic classiher for the mul¬ 
tirate system closely follows the development for the single-channel and multichannel 
cases developed in [Ref. 22] and [Ref. 21]. It is necessary, however, to make modih- 
cations to some of the parameters, to account for the periodic nature of the multirate 
system. Thus, some parameters that are constant in the single-channel and multi¬ 
channel cases become periodically time-varying variables in the multirate case. 

To start the development, assume that the signals associated with both classes 
are jointly Gaussian and that the observations for all channels are collected and 
ordered such that x„ contains all observations from time 0 to n. In addition, the 
observation vector (5.1) can be written recursively as 

(5.20) 


x„ = 


^n—1 

(k) 
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where the superscript k is used to indicate that the size of the observation vector 
is a periodically time-varying parameter since the number of observed channels 
at any time step varies periodically. If the size of (the number of observations at 
time n) is dehned as and the total number of observations from time 0 to n — 1 
is dehned as then the total number of observations at time n can be expressed 

as 

= (5.21) 

As an example, consider the two-channel system where the hrst channel is 
decimated by a factor of 2 and the second channel is decimated by a factor of 3, as 
shown in Fig. 5.1. The system period as dehned in Chapter II is iF = 6. Thus, since 
n = k (mod K), the system phase, fc, varies periodically between 0 and 5. 


n 

X, 



29 30 31 32 33 34 35 36 37 38 


X 


2 


o 


o 


X 


29 



o 


Figure 5.1. Two-channel Multirate System. 


This can be illustrated more specihcally as follows (refer to Fig. 5.1). Assum¬ 
ing that observations through n = 29 have been collected, the associated observation 
vectors for time steps n = 30, 31,..., 35 would be denoted by 



X 29 


X 30 


X 31 

X 30 = 

^( 0 ) 

X 31 = 

^( 1 ) 

X 32 = 



2i30 


2i31 


±32 


X 32 


X 33 


X 34 

X 33 = 

^(3) 

X 34 = 

^(4) 

X 35 = 



2i33 


2i34 


2i35 
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where the vectors of observations at these times are 


^(0) _ 
±30 — 


:ri[30] 

a;2[30] 


^ 31 ^ — empty ^2 — a^i[32] 


^ 33 ^ = 2^2 [33] ^ 34 ^=a;i[ 34 ] ^ 5 = empty 

The corresponding number of observations, , at these times are 


c(*^) 9 cW fi c'(2) 1 

^30 ~ ^ ^31 ~ ^32 — 


c(3) _ 1 c(4) _ 1 c(5) _ ri 

^33 ~ ^34 — ^35 ~ 


For any given time ‘n’, the probability density function for the collection of 
observations is given by the multirate Gaussian form 


/x„(Xn) — 


(27r)^»/2|C„|V2 


g-i(x„-m„)^C„bxn-m„) 


(5.22) 


where m„ = T{x„} is the mean vector and C„ = £^{(x„ — m„)(x„ — m„)^} is the 
covariance matrix for the observations. The mean vector and covariance matrix can 


then be written as 




(5.23) 


C„ = 


TD 

'^n—1 

■p(fc)^ Yi(fc) 

iv,n 


(5.24) 


where the partitioning corresponds to the partitioning dehned in (5.20). By using a 
formula for matrix partitioning [Ref. 69], the inverse covariance matrix, C“^, can be 
represented in terms of the inverse covariance matrix, as 


c-i = 


I c-\ 0 


I 0 


0 ^ I 


0 ^ I 


(5.25) 
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or equivalently as 


o' 

+ 

(fc) 

n 

I 

(5.26) 

0^ 0 


I 





where and are defined as 



C 


-1 
n—1 






( 5 . 27 ) 

(5.28) 


Since the determinants of the upper and lower triangular matrices of (5.25) are unity, 
the determinant of C“^ is given by 


|C;^I = IC 


-1 
n—1 




or 



(5.29) 


The partitionings of (5.1) and (5.23) through (5.29) provide all the necessary relations 
for the recursive computation of the density function (5.22) at each time step. 

Now that the probability density function for the multirate system has been de¬ 
scribed, the probability density function for a given class within the multirate system 
can be discussed. Consider observations from a known class i with a corresponding 
mean vector 


m 


(b 


m 


m 


(b 

n—1 

{i){k) 


(5.30) 


and covariance matrix 



p((b 

^n-l 

■p (O(^) 

-Tvn 

■p (b(^)^ 



(5.31) 


where the superscript {i) denotes association with class i and the partitioning corre¬ 
sponds to the partitioning defined in (5.23) and (5.24). The vector of observations 
with the mean removed can be denoted as 


xh) — X 


— 


(5.32) 


92 



or 


xi*) = 


xh) 


Aim 



(*) 

Xn-l - 
(fc) _ (*)(fc) 

vOn jJJ±n 


(5.33) 


By using the form (5.22) and substituting the relations (5.21), (5.32), (5.33), (5.26) 
and (5.29), the probability density function for class i at time n can be written as 


/x„|_H'i(Xn|hrj) — 




(27r)(^»)/2|ci'V/' 


1 f 


X 


(27r)^«-i/2|C^l,|V2 (27r)^»''V2|E«('=)|i/2 


X e 


The hrst term on the right represents /x„_i|Hi while the second term on the right 
represents a conditional probability density function /x„|x„_iHi- In other words, (5.34) 
can be written as 

/x„|Hi = /x„_i|_ffi ■ /x„|x„_iHi i = 1,2 (5.35) 

where 

= T—T-(5-36) 

and 

1 






(27r)5»''V2|Ei*)(^V/^ 


(27r)^»-i/2|C^li|V2 

(Jim \ 

g 2 V—" ^n-1 j y^ri '-^n X„_j J 


Using these results, the statistic h„(x„) dehned in (5.10) becomes 

h„u„) =-2hi I ■ i--'"-'"' 


_ 2 In ^ 1^1 (x^ 11 ) 2 In |xn—i 

/x„_i|Ho(^«-l l-^o) /x„|x„_iHo 


(5.37) 


(5.38) 
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This can be written as 


where 



(5.39) 

^In 

fxJXn-lHo 

(5.40) 


Snbstituting (5.37) into (5.40) resnlts in the expression 


1 4 - 

p(fc)(l)' 

" (1) 42 (AO) 
^n-l) {An 

(fc)(o)- 

^n-l) 1 


|Xn—1 ) — 

\^n 



1 111 - 



(5.41) 


Since the hrst observation vector is labeled xq and the statistic ho(xo) = 
h_i(x_i) + Aho(x® |x_i), the SPRT is hrst initialized with X-i = 0 and the statistic 
h_i(x_i) is empty. Thereafter, is compnted from (5.39) and (5.41) recursively. 
This statistic is used in the decision rule (5.9) to perform the classihcation. 


C. CLASSIFICATION USING LINEAR PREDICTION 


The equations for the multirate sequential classiher were developed assuming 
that the observations were jointly Gaussian. The sequential classiher, however, can 
also be interpreted in terms of linear prediction. By rearranging (5.27) and (5.28) 
into 

(5.42) 

Ef = E<,‘'-Rf"GW, (5.43) 


these two equations can be combined in matrix form to write 


y-i(fc) 

i 

£g 


I 


1 

H 

_1 

1 

1 

u 


1 

$e 

0 

1 


-1 

o 


(5.44) 


These equations bear a striking similarity to the Normal equations of linear 
prediction (see (4.2)). Indeed, the matrix —Gn'^ can be interpreted as the prediction 
hlter coefficients for predicting the new observations, and is the prediction error 
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covariance matrix. More explicitly —and relate to the parameters in (4.9) 
by 


and 


Aq 

Ai 

Aq 


E = 


I 

1 0 1 . 



1 I 1 • 



1 1 . 


-L'n 1 

X • • • 

X 

0 1 

iLn \ • • • 

X 

; 1 

: 1 

X 

0 1 

0 1 ... 

1 E^°^ 

1 -L'n 


(5.45) 


(5.46) 


The block form in (5.45) and (5.46) may be a little misleading. In most cases, 
the width of the blocks is just one column and the becomes a “1.” The matrices 
in (5.45) and (5.46) have been represented in block form to allow for possible multiple 
observations at each time step. The parameters —and En^^ thus represent the 
periodically time-varying parameters used to predict each new observation (or set of 
observations) within a block. 

When a sufficiently high order Q is chosen for the prediction, the parameters 
^ and E(f ^ do not depend explicitly on the time variable “n” but only on the phase 
as indicated by the superscript k. 


D. ALGORITHM SUMMARY 

This section summarizes the multirate sequential classiher algorithm for ease 
of implementation. The recursive classihcation algorithm can be divided into two 
parts: 1) training and 2) testing. 

In the hrst part, the system is trained to distinguish between two classes using 
a priori information. Using known class training data, the classiher parameters (mean 
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vectors, prediction coefficient matrices and prediction error variance matrices) for each 
class are estimated. In addition, the thresholds used in classihcation are determined 
and adjusted experimentally. These parameters are summarized in Table 5.1. 


Table 5.1. Summary of Classiher Training Parameters 


Class i Parameters: 

mean vectors, = Si[xn '^] 

i = 0,l 

Prediction coefficients, 


Prediction Error Variances, 

Classihcation Thresholds: 

A (Class 1) and B (Class 0) 


In the second part new (test) data is presented to the classiher. The associated 
class errors and ) for the given signal, x„, are calculated and used 

along with the class prediction error variances to calculate the SPRT statistic, /i„, 
for comparison with the class thresholds. When the SPRT statistic crosses a class 
threshold, the classiher declares the classihcation; otherwise, the process is repeated. 
These steps are shown in Table 5.2. 


Table 5.2. Steps in Sequential Classihcation 
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The performance of the classiher can be evaluated with two key measures of 
performance: the number of correct classihcations and the length of data (length of 
time) required to make a classihcation. By varying different properties of the classiher 
and environment (e.g., signal-to-noise ratio), the effects of these different properties on 
the classiher can be evaluated using the two performance measures. The remainder 
of this chapter describes the simulations and results for varying conditions of the 
classiher and environment. 

E. SIMULATION SETUP AND PARAMETERS 

In order to test the proof of concept for the multirate sequential classiher 
multiple simulations were conducted. The test data consisted of recordings of sound 
from one propellor plane and three diherent A-10 jet aircraft, with noise added. The 
parameters that were varied during the simulations were the signal-to-noise ratio 
(SNR) associated with the training and target data sequences, and the length of 
the training sequences used to develop the classiher parameters. Simulations were 
conducted to compare the single-channel, multichannel and multirate cases. 

1. Test Data Characteristics 

This subsection describes the characteristics of the test data. The original 
audio data sequences were recorded at 44,100 Hz. The data was initially low-pass 
hltered and downsampled by a factor of 10 to achieve a sampling rate of 4,410 Hz. 
(The observed data for the simulation is then further decimated.) A spectral estimate 
was then computed for each of the four data sequences prior to adding any additional 
noise, in order to determine any unique characteristics. Figure 5.2 shows the spectral 
plots of each sequence and Table 5.3 summarizes the dominant frequencies, and a 
secondary frequency if noticeable, that were measured from the spectral plots. 

From Fig. 5.2, it can be seen that the propellor plane has secondary frequency 
at a lower frequency than the peak frequency of 192.0 Hz. A secondary frequency is 
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Figure 5.2. Spectral Plots 


Table 5.3. Spectral Information for Data Sequences 



Propellor Plane 

A-10 Jet No. 1 

A-10 Jet No. 2 

A-10 Jet No. 3 

Dominant 
Freq (Hz) 

192.0 

145.0 

82.5 

84.6 

Secondary 
Freq (Hz) 

106.9 

N/A 

N/A 

210.4 
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also seen on the A-10 Jet No. 3 at 210.4 Hz. The A-10 Jet No. 1 has many peaks close 
to each other. However, the magnitude of these peaks diminish as the distance from 
the main peak increases. Thus they have not been considered as secondary peaks. 

In addition to the spectral estimate of the test data, the autocorrelation of each 
data sequence was calculated to examine any signihcant differences in correlation. 
The autocorrelation plots were used to assist in selecting appropriate hlter orders for 
calculating the prediction orders used in the classiher. One can see from Fig. 5.3 that 
the three A-10 Jet planes have similar, but not identical, correlation structure within 
the hrst ten samples (dashed lines in Fig. 5.3.) This similar structure led to choosing 
a prediction hlter order larger than ten, as discussed in the next subsection. 




Lag Lag 

(a) Autocorrelation for Propeller Plane (b) Autocorrelation for A-10 Jet No. 1 




Lag Lag 

(c) Autocorrelation for A-10 Jet No. 2 (d) Autocorrelation for A-10 Jet No. 3 


Figure 5.3. Estimated Autocorrelation Function for Aircraft Data at 4,410 Hz sam¬ 
pling rate. Dashed lines show ±10 lag values. 
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2. Simulation Models 

For the simulations conducted, three scenarios were considered. These scenar¬ 
ios considered classihcation using a single-channel system, the multichannel system 
and the multirate system. All scenarios began with the recorded data sampled at 
4410 Hz. For the single-channel scenario, the recorded data was decimated by a fac¬ 
tor of two, resulting in a sampling rate of 2205 Hz. This was done to ensure that 
the highest rate in all three scenarios would be the same. Additive white Gaussian 
noise with a mean of zero and a variance of cr^ was then added to achieve a desired 
Signal-to-Noise ratio set by the simulation parameters. The proper selection of the 
value of a'^ is described below. The single-channel scenario is shown in Fig. 5.4. 



Channel 1 
(2205 Hz) 


Figure 5.4. Observation Model for the Single Channel Scenario. 


The simulation model for the multichannel scenario is similar to the single¬ 
channel scenario, with the exception that after the recorded data is decimated by a 
factor of two, the signal is split into two channels. Noise was then added to each 
channel to achieve a desired SNR. The model for the multichannel scenario is shown 
in Fig. 5.5. 

For the multirate scenario. Fig. 5.6, the high rate channel was decimated by a 
factor of 2 and the low rate channel was decimated by a factor of 3. Noise was then 
added to each channel to achieve a desired signal-to-noise ratio. 

In order to produce the desired SNR for these models, the noise power must 
be calculated using the desired SNR in decibels and the estimated power of the 
appropriate channel. If the output of the decimator is s[mx], where rrix is the time 
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Noise 1 



Channel 1 
(2205 Hz) 


Channel 2 
(2205 Hz) 


Noise 2 


Figure 5.5. Observation Model for the Multichannel Scenario. 


Noise 1 



Channel 1 
(2205 Hz) 


Channel 2 
(1470 Hz) 


Noise 2 


Figure 5.6. Observation Model for the Multirate Scenario. 
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index for the decimated signal, f][mx] is the noise signal, and the observed signal x[mx] 

is 

x[mx] = s[mx] +r]['fT^x], 

then the noise power, can be calculated as 

rrix-l 

^ E 

Cr2 = __ 

v ]^g-SNR/10 

Once the noise power has been calculated, this value can be used to generate the 
appropriate noise signal, ri[mx]. 

3. Prediction Coefficients 

When calculating the prediction coefficients for the training data, an appro¬ 
priate hlter order has to be chosen. In general, the hlter order should be large enough 
to exploit differences in correlation between the two classes of signals being tested. 
It should be small enough, however, so that it does not add excessive computational 
burden to the simulation and does not become sensitive to small artifacts of the train¬ 
ing data. In addition, as the hlter order increases, the determinant of the prediction 
error covariance, which measures the quality of the prediction, typically approaches 
a value such that any further increases in hlter order do not result in an appreciable 
improvement in performance. 

From Fig. 5.3, it can be seen that the propellor plane is uniquely diherent in its 
correlation structure, but the three A-10 jets have similar structures for lags smaller 
then ten. Thus choosing an order for the linear prediction hlter less than ten would 
probably not provide enough uniqueness in the three A-10 jet parameters to be able 
to adequately classify the diherent planes. However, the correlation shapes at lags 
larger than ten are sufficiently dissimilar to exploit these diherences in correlations. 
Although larger hlter orders could have been used, the hlter order was chosen to be 
40 to minimize excess computations. 


102 



The prediction coefficients and prediction error variances were calculated in the 
manner described in Chapter IV, using the multichannel Levinson recursion. Since 
the decimation rates for the two channels in the multirate scenario are two and three, 
the system period, K, is equal to 6. In addition, there are a total of hve observations 
that occur within one system period, three for the high-rate channel and two for the 
low-rate channel. The sampling pattern is shown in Fig. 5.7. It is assumed the data 
is aligned so that at the initial observation, both channels produce observations. 


System Period K = 6 

/ -\ 


0 12345678 


Channel 1 


Channel 2 


K,=2 


O 


_ ) -^ 

if 2 = 3 time 


Figure 5.7. Illustration of Decimation for the Multirate Scenario. 


4. Training Data 

When calculating the prediction coefficients, a set of training data for each 
class had to be used. This data consisted of a portion of the total data signal used 
for each class, as shown shaded in Fig. 5.8. The length of this portion was varied 
according to the scenario parameters. The rest of the signal was then used for the 
testing portion of the simulation. 

During the testing phase, the simulation starts at the beginning of the test 
portion of the input signal and continues to take observations until a classihcation is 
made. This is considered to be one “trial.” A second trial is then begun starting with 
the observation corresponding to the data point just after the data point at which 
the hrst classihcation was made, as shown in Fig. 5.8. This process is continued 
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until 50 trials are conducted. The length of the data was sufficient to ensure that the 
simulation was able to conduct 50 trials before reaching the end of the data. 



Start of 
Classifier Data 

Figure 5.8. Training Data. 


Although the sequential classiher can in theory classify data based on as little 
as one single observation, in practice it is better to allow some minimum number of 
observations to be collected before a classihcation is attempted. In these experiments, 
this minimum data length was linked to the length of the linear prediction filter. In 
particular, the algorithm was designed not to classify the target until the number of 
observations exceeded the hlter order by one observation. It was decided to implement 
the algorithm this way to avoid any errant classifications that might occur while the 
algorithm collects the observations needed to make the predictions at full order. 


F. SIMULATION RESULTS AND ANALYSIS 

The results for the multirate scenarios are presented in the following subsec¬ 
tion; the key parameters are analyzed in the last subsection. 
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1 . 


Multirate Scenario Results 


Since there were four classes of data (one propellor and three jet planes) the 
classiher was tested using two classes at a time. This resulted in six sets of experi¬ 
ments, whose results are listed in Tables 5.4 through 5.9 of this section. 

In order to make comparisons, specihc parameters were varied for the simula¬ 
tions. These were the signal-to-noise ratio (SNR) of the target data and the length 
of data sequences used to calculate class parameters. In all cases, the training data 
was assumed to be noiseless but white noise was added to the test data to achieve a 
desired SNR. 

Table 5.4. Classihcation Results: Propellor vs. A-10 Jet No. 1 


SNR 

Data 

Training 

Length 

Propellor 

No. of Correct 

Classifications 
(out of 50) 

Jet No. 1 

No. of Correct 

Classifications 
(out of 50) 

Average Time 
to Classify 
Propellor 
(msec) 

Average Time 
to Classify 

Jet No. 1 
(msec) 

40 dB 

6250 

50 

50 

9.52 

9.52 

2500 

50 

50 

9.52 

9.52 

20 dB 

6250 

50 

50 

9.52 

9.59 

2500 

50 

39 

9.52 

9.52 

10 dB 

6250 

50 

21 

9.52 

9.52 

2500 

50 

10 

9.52 

9.52 


Table 5.4 shows the classihcation results for the propellor plane versus the 
A-10 Jet No. 1. This table lists both the number of correct classihcations for each 
target as well as the average time to classify the target. The average time of 9.52 msec 
listed in the table corresponds to the minimum observation condition of 41 samples 
(see Subsection V.E.4). It can be seen from Table 5.4 that the classiher makes its 
decision almost immediately after reaching the full hlter order under all cases. Thus, 
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for the case of the propellor plane versus Jet No. 1, the advantage of a sequential 
classiher is not realized. This is not true for all the test scenarios. 

From Table 5.4, one can see that the algorithm is able to classify the propellor 
under all conditions, and almost immediately after collecting enough data to reach 
full order. In an environment with an SNR of 40 dB, the classiher is able to classify 
the propellor plane nearly perfectly. As the SNR is reduced to 20 dB, the success 
rate reduces to a little below 80% for the short training data length and under 15% 
for the high training data length. With the middle length training data classihcation 
is perfect. Proportionately similar reductions occur as the SNR is dropped further to 
10 dB. The results degrade signihcantly at an SNR of 10 dB. Namely, the classiher 
tends to classify both signals as a propellor plane; this puts into question whether 
the classiher appropriately classihes the propellor plane at 10 dB or defaults to the 
propellor plane because of the signal degradation. In almost all cases classihcation 
occurred with the minimum allowable number of observations. For SNR below 5 
dB, classihcation degraded severely. Therefore, the results are not presented in these 
tables. 


Table 5.5. Classihcation Results: Propellor vs. A-10 Jet No. 2 


SNR 

Data 

Training 

Length 

Propellor 

No. of Correct 

Classifications 
(out of 50) 

A-10 No. 2 

No. of Correct 

Classifications 
(out of 50) 

Average Time 
to Classify 
Propellor 
(msec) 

Average Time 
to Classify 
A-10 No. 2 
(msec) 

40 dB 

6250 

50 

50 

9.56 

9.52 

2500 

50 

50 

9.53 

9.52 

20 dB 

6250 

50 

36 

9.52 

9.53 

2500 

50 

8 

9.52 

9.53 

10 dB 

6250 

50 

2 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 


The results of classihcation for the propellor plane and the A-10 Jet No. 2 
are contained in Table 5.5. Similar to the testing between the propellor plane and 
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the A-10 Jet No. 1, the algorithm was able to correctly classify the propellor plane 
nearly perfectly under all conditions simulated. In addition, the classihcation occurred 
almost immediately after reaching the minimum number of observations allowed. The 
number of correct classihcations for Jet No. 2 was slightly less than that for Jet No. 1. 
As with Jet No. 1, the performance of Jet No. 2 dropped appreciably as the SNR 
was reduced to 20 dB. Under the conditions with an SNR of 10 dB, the classiher 
performed poorly when trying to classify Jet No. 2. 


Table 5.6. Classihcation Results: Propellor vs. A-10 Jet No. 3 


SNR 

Data 

Training 

Length 

Propellor 

No. of Correct 

Classifications 
(out of 50) 

A-10 No. 3 

No. of Correct 

Classifications 
(out of 50) 

Average Time 
to Classify 
Propellor 
(msec) 

Average Time 
to Classify 
A-10 No. 3 
(msec) 

40 dB 

6250 

50 

50 

9.52 

9.52 

2500 

50 

50 

9.52 

9.52 

20 dB 

6250 

50 

45 

9.52 

9.53 

2500 

50 

35 

9.52 

9.53 

10 dB 

6250 

50 

2 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 


For the experiments conducted between the propellor plane and the A-10 Jet 
No. 3 shown in Table 5.6, the algorithm once again was able to classify the propellor 
plane one hundred percent of the time and with the minimum number of observations 
allowed. The results for Jet No. 3 were somewhat between the result for the other 
two A-10 jets. For the simulations with an SNR of 40 dB, the algorithm was able 
to classify Jet No. 3 perfectly. As the SNR was reduced to 20 dB, the performance 
dropped to between 70% and 90%. The performance then signihcantly dropped to 
almost no correct classihcations with an SNR of 10 dB. 

Besides comparing the jet planes to the propellor plane, it was also of interest 
to see if the classiher could discriminate between the diherent types of jet planes. 
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Table 5.7. Classification Results: A-10 Jet No. 1 vs. A-10 Jet No. 2 


SNR 

Data 

Training 

Length 

A-10 No. 1 

No. of Correct 

Classifications 
(out of 50) 

A-10 No. 2 

No. of Correct 

Classifications 
(out of 50) 

Average Time 
to Classify 
A-10 No. 1 
(msec) 

Average Time 
to Classify 
A-10 No. 2 
(msec) 

40 dB 

6250 

49 

47 

10.01 

9.94 

2500 

11 

50 

9.68 

9.52 

20 dB 

6250 

43 

23 

9.68 

9.67 

2500 

19 

42 

9.65 

9.54 

10 dB 

6250 

43 

16 

9.54 

9.57 

2500 

15 

37 

9.52 

9.52 


This is obviously a more difficult problem than discriminating between a jet plane 
and a propellor plane. 

When comparing the different jets, the results were more balanced than those 
seen comparing the jets with the propellor plane. For the experiments conducted 
comparing Jet No. 1 and Jet No. 2, shown in Table 5.7, percentages of correct clas¬ 
sification greater than eighty percent were typical for about half of the scenarios. 
The rest of the results ranged from 0% to about 50%. The time needed for the al¬ 
gorithm to make a classification decision is somewhat higher than the time needed 
when comparing a jet plane to the propellor plane. 

Table 5.8 contains the results for the experiments conducted comparing the 
A-10 Jet No. 1 with Jet No. 3. From this table one can see that the algorithm was 
able to successfully classify Jet No. 1 correctly for the environments with SNRs of 20 
dB and 10 dB, while Jet No. 3 performed poorly under these conditions. The results 
were opposite for the case of an SNR of 40 dB. Classifications in the 20 dB and 10 dB 
cases were made almost immediately after reaching the minimum allowable length. 
For the 40 dB case the classification time was significantly longer. 

Results for comparing the A-10 Jet No. 2 and the A-10 Jet No. 3 are contained 
in Table 5.9. From this table one can see that nearly perfect results for classifying 
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Table 5.8. Classification Results: A-10 Jet No. 1 vs. A-10 Jet No. 3 


SNR 

Data 

Training 

Length 

A-10 No. 1 

No. of Correct 

Classifications 
(out of 50) 

A-10 No. 3 

No. of Correct 

Classifications 
(out of 50) 

Average Time 
to Classify 
A-10 No. 1 
(msec) 

Average Time 
to Classify 
A-10 No. 3 
(msec) 

40 dB 

6250 

11 

50 

10.56 

9.52 

2500 

5 

50 

9.53 

9.52 

20 dB 

6250 

50 

1 

9.52 

9.66 

2500 

50 

2 

9.52 

9.52 

10 dB 

6250 

50 

0 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 


Table 5.9. 


Classification Results: A-10 Jet No. 2 vs. 


A-10 Jet No. 3 


SNR 

Data 

Training 

Length 

A-10 No. 2 

No. of Correct 

Classifications 
(out of 50) 

A-10 No. 3 

No. of Correct 

Classifications 
(out of 50) 

Average Time 
to Classify 
A-10 No. 2 
(msec) 

Average Time 
to Classify 
A-10 No. 3 
(msec) 

40 dB 

6250 

40 

49 

12.48 

11.40 

2500 

50 

49 

13.39 

11.65 

20 dB 

6250 

49 

4 

9.53 

9.56 

2500 

50 

0 

9.52 

9.52 

10 dB 

6250 

50 

0 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 
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Jet No. 2 were obtained. Classifications occurred almost immediately after reaching 
full order for the 20 dB and 10 dB cases. While more time was needed for the 40 
dB cases. For Jet No. 3 nearly perfect classification results were obtained for the 40 
dB cases. This performance fell to nearly 0% for the 20 dB ad 10 dB cases. Similar 
to the experiments with the propellor plane and Jet No. 1, as the SNR is dropped 
significantly, degradation only in Jet No. 3 tends to show that the classifier is skewed 
towards Jet No. 2 questioning whether the results for Jet No. 2 at SNRs of 20 dB 
and 10 dB are reliable. 

In another set of experiments, the classifier designed using the propellor and 
Jet No. 1 was tested on data from all four classes (propellor. Jet No. 1, Jet No. 2 and 
Jet No. 3). When tested on Jet No. 2 and Jet No. 3, the classification was considered 
to be correct if the classifier chose the “Jet” category, and incorrect if the classifier 
classified these jet planes as propellor planes. 

Results of these experiments are listed in Table 5.10. From these results, it 
can be seen that the classifier was able to correctly classify the A-10 Jets No. 2 and 
No. 3 as jets perfectly for the 40 dB scenarios. When the SNR was reduced to 20 
dB, the A-10 Jet No. 3 was classified nearly perfectly, while the A-10 Jet No. 2 had 
at most a classification rate of about 60%. Under the 10 dB conditions, the classifier 
performed poorly when trying to classify the A-10 Jets No. 2 and No. 3. 

2. Analysis of Single-channel vs. Mnltirate vs. Mnltichannel 
Simulations 

Experiments identical to those conducted for the multirate data were con¬ 
ducted for the single channel and multichannel scenarios. In most situations, it was 
possible to adequately classify the target in these other scenarios. Therefore, to com¬ 
pare the effectiveness of the methods, the average number of data samples needed to 
make a classification were also compared. Tables 5.11 through 5.13 show results for 
some of the experiments for an environment with an SNR of 40 dB. Table 5.11 shows 
that under these conditions, the classifier was able to correctly classify the propellor 
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Table 5.10. Classification Results: Class 0 = Propeller and Class 1 = A-10 Jet No. 1 


SNR 

Data 

Training 

Length 

Propellor 

No. of Correct 

Classifications 
(out of 50) 

Jet No. 1 

No. of Correct 

Classifications 
(out of 50) 

Jet No. 2 

No. of Correct 

Classifications 
(out of 50) 

Jet No. 3 

No. of Correct 

Classifications 
(out of 50) 

Average Time 
to Classify 
Propellor 
(msec) 

Average Time 
to Classify 

Jet No. 1 
(msec) 

Average Time 
to Classify 

Jet No. 2 
(msec) 

Average Time 
to Classify 

Jet No. 3 
(msec) 

40 dB 

6250 

50 

50 

50 

50 

9.52 

9.52 

9.52 

9.52 


2500 

50 

50 

30 

50 

9.52 

9.52 

10.24 

9.52 

20 dB 

6250 

50 

50 

32 

47 

9.52 

9.52 

9.57 

9.53 


2500 

50 

39 

3 

41 

9.52 

9.75 

9.55 

9.54 

10 dB 

6250 

50 

21 

2 

9 

9.52 

9.52 

9.52 

9.52 


2500 

50 

10 

0 

10 

9.52 

9.52 

9.52 

9.53 




and Jet No. 2 in all these scenarios. In addition, it can be seen that the time needed 
to make the classihcation is roughly the same for the single-channel, multichannel 
and multirate cases. 

Table 5.11. Classihcation Results for Propellor and Jet No. 2 


Scenario 

Data 

Training 

Length 

Propellor 

No. of Correct 
Classihcations 
(out of 50) 

Jet No. 2 

No. of Correct 

Classihcations 
(out of 50) 

Average Time 
to Classify 
Propellor 
(msec) 

Average Time 
to Classify 

Jet No. 2 
(msec) 

Single-Channel 

6250 

50 

50 

10.12 

9.56 


2500 

50 

50 

9.56 

9.52 

Multichannel 

6250 

50 

50 

9.52 

9.52 


2500 

50 

50 

9.76 

10.28 

Multirate 

6250 

50 

50 

9.53 

9.52 


2500 

50 

50 

9.52 

9.53 


In table 5.12, it can be seen that ability to distinguish between Jet No. 1 and 
Jet No. 2 becomes more difficult than between the propellor plane and jets for all 
three scenarios. The multirate classiher, however, is able to make its classihcations 
in a shorter amount of time. 

In the example of the A-10 Jet No. 2 versus the A-10 Jet No. 3, Table 5.13, 
the classihers again performed fairly well. Although, the amount of time necessary 
to make classihcation decisions for all three cases was slightly higher than for some 
of the other simulations, the multirate classiher on average needed less time to make 
its classihcation. 

From these tables it can be seen that the multirate case has some advantage 
over the single-channel case. In low noise environments where the single-channel, 
multirate and multichannel scenarios all had a high success rate of proper classihca¬ 
tion, the multirate scenario was able to classify the proper signal in a shorter amount 
of time. It was noticed that the multichannel channel case at times needed more time 
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Table 5.12. Classification Results for Jet No. 1 and Jet No. 2 



Data 

Jet No. 1 

Jet No. 2 

Average Time 

Average Time 

Scenario 

Training 

No. of Correct 

No. of Correct 

to Classify 

to Classify 


Length 

Classihcations 

Classifications 

Jet No. 1 

Jet No. 2 



(out of 50) 

(out of 50) 

(msec) 

(msec) 

Single-Channel 

6250 

50 

36 

9.90 

35.78 


2500 

45 

49 

15.96 

22.78 

Multichannel 

6250 

50 

35 

9.89 

32.42 


2500 

46 

50 

15.37 

22.51 

Multirate 

6250 

49 

47 

10.02 

9.94 


2500 

11 

50 

9.68 

9.52 


Table 5.13. Classification Results for Jet No. 2 and Jet No. 3 



Data 

Jet No. 2 

Jet No. 3 

Average Time 

Average Time 

Scenario 

Training 

No. of Correct 

No. of Correct 

to Classify 

to Classify 


Length 

Classihcations 

Classihcations 

Jet No. 2 

Jet No. 3 



(out of 50) 

(out of 50) 

(msec) 

(msec) 

Single-Channel 

6250 

35 

50 

20.10 

10.74 


2500 

50 

49 

15.13 

13.08 

Multichannel 

6250 

37 

50 

17.29 

11.31 


2500 

50 

48 

12.05 

15.18 

Multirate 

6250 

40 

49 

12.48 

11.40 


2500 

50 

49 

13.48 

11.65 
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to make a classification decision than the multirate case. Although there is more 
data available in the multichannel scenario, the offset observations of the multirate 
scenario adds more statistical information to the solution. 

3. Key Parameters for the Multirate Case 

There were three main characteristics of the data that were analyzed for their 
effect on the classihcation performance. They were the dominant frequency of the 
target data, the signal-to-noise ratio and the training data length. In addition, the 
performance between single-channel, multirate and multichannel observations were 
compared. Each is described individually in the following subsubsections, however, 
their impact is not independent of each other. So, if noticeable, their relative impact 
on each other is discussed. 

a. Dominant Frequency 

While analyzing the data, it became apparent that the classiher per¬ 
formed much better when one target data class had a dominant frequency higher than 
that of the other class to which it was being compared. This was especially noticeable 
in the high noise environment scenarios. Under the high noise conditions, the jets 
were not classihed correctly with any degree of success against the propellor plane. 
Only under the low noise condition [SNR = 10 dB) did the jets display moderate 
success with classihcation against the propellor plane. Likewise, Jet No. 1 had a 
higher correct classihcation rate when attempting to classify against Jets No. 2 and 
No. 3. And, Jet No. 2 performed better against Jet No. 3. 

b. Signal-to-Noise Ratio 

By analyzing the simulation data, it became noticeable that under low 
noise conditions {SNR = 40 dB), the three A-10 jets were classihed properly against 
the propellor plane. As the training sample length was reduced, the performance 
degraded dramatically. In addition, the higher noise conditions reduced the 
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performance of the jets significantly when compared to the propellor plane. It appears 
that the dominant frequency has a large effect on the outcome of the classiher, in any 
condition other than ideal. 

c. Training Data Length 

As might be anticipated, the longer the training data length, the better 
the classiher typically performed. Not only did the classihcation rate generally in¬ 
crease, but the number of samples needed to make a proper classihcation was reduced. 
For the audio data used, however, when the size of the training data was increased to 
12500 samples the performance of the classiher actually degraded. It was determined 
that at this length the data available for testing became statistically diherent enough 
from the training data and caused a degradation in the performance of the classihers. 

G. CONCLUSION 

In this chapter, it was shown that with some modihcations to the single¬ 
channel classiher, one can implement a sequential classiher that accepts multirate 
data. In addition, using the multirate Levinson recursion developed in Chapter IV; 
one can derive the needed class parameters for the multirate sequential classiher. The 
multirate sequential classiher can handle multiple channels of observations occurring 
at diherent sampling intervals. 

Simulations were conducted to analyze the ehectiveness of the multirate se¬ 
quential classiher compared with a single-channel classiher and a multichannel clas¬ 
siher where all channels were sampled at the same rate. It was seen that under low 
noise environments all three classihers performed well, but the multirate classiher was 
able to make classihcation decisions in a shorter amount of time due to the added 
information from the ohset observations. The amount of addiitonal information con¬ 
tained in the multirate classiher is dependent upon the disparity between channel 
sampling rates. 
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VI. CONCLUSIONS AND FUTURE WORK 


The goal of this work was to further develop the theory and applications of 
statistical multirate signal processing. A summary of the work and contributions of 
this thesis is provided here. This is followed by a section of suggestions for further 
research. 

A. CONCLUSIONS 

In Chapter II, the foundation necessary to describe multirate systems is es¬ 
tablished. Part of this foundation requires dehning some key terms applicable to 
multirate systems, such as fundamental rate, system period and system phase. These 
terms are used to explicitly describe the periodic nature of the multirate system 
and allow for the appropriate selection of values for the periodic components in the 
multirate system. The dehnitions for describing the statistical characterizations of 
multirate signals and building blocks (e.g., decimator and expander) are provided. 
Compact matrix forms of the decimator, expander and hlters, which take advantage 
of Kronecker product relations, are presented. Finally, the LPTV hlter bank of [Ref. 
17] is generalized to apply to systems with different input and output rates. 

In Chapter III, the explicit direct form of the multirate optimal estimator is 
developed. This optimal estimator uses multiple-input signals observed at different 
sampling rates to estimate a desired signal. In addition, a recursive innovations 
form of this multirate optimal estimator is derived. This innovations form of the 
optimal estimator separates the direct form optimal hlters into modihed optimal 
hlters and cross hlters. These cross hlters remove any information from one signal 
that is contained in the other signals, in an ordered fashion. In essence, a new set 
of input signals that are mutnally orthogonal are derived and used as the inputs to 
the modihed optimal hlters. Using the recursive form of the optimal hlter allows 
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one to calculate the relative change in performance in the optimal estimator when 
input signals are added or removed from the system. Through simulations, it is 
demonstrated that optimal hltering of multiple multirate channels for an AR process 
and a sinusoidal process can provide improved performance over optimal hltering 
using a single channel, even if the secondary channel when used by itself has high 
error variances. 

In Chapter IV, the optimal hltering problem is specialized for the case of 
optimal multirate linear prediction. The multirate Normal equations are derived 
for a multirate system with multiple input and output signals, which are observed at 
diherent sampling rates. For a multirate system with a system period K, there are up 
to K distinct sets of prediction coefficients and error covariance matrices that apply in 
a periodic fashion. An efficient method for calculating the multirate linear prediction 
coefficients and error variances is developed through the use of the multichannel 
Levinson recursion and generalized triangular UL factorization. A specihc algorithm 
for the generalized triangular factorization is included in Appendix E. 

In Chapter V, a multirate sequential classiher is derived starting from the 
basic theory of sequential hypothesis testing. It is shown that classiher parameters 
needed for implementing the multirate sequential classiher are the same as those for 
multirate linear prediction. A multirate sequential classiher is then implemented and 
tested using audio hies of a propellor plane and three A-10 jet aircraft. The experi¬ 
ments tested the classiher performance in selecting between the propellor plane and 
jet aircraft as certain system parameters are changed. These parameters consisted of 
the signal-to-noise ratios of the observed signal and the length of the training data. 
In addition, the performance of the multirate classiher is compared to that of the 
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single-channel and multichannel classihers using similar data. These experiments 
showed that improved performance can be obtained by using two channels of obser¬ 
vation data even if one channel is sampled at a lower rate. 

B. TOPICS FOR FURTHER INVESTIGATION 

Multirate signal processing continues to be an active area of research and many 
opportunities exist for further research, particularly in the area of statistical signal 
processing. The optimal multirate hlter developed in this dissertation is an FIR 
hlter. The inhnite impulse response (HR) hlter was not considered here; however, 
multirate theory applies equally well to the HR case and also would provide signihcant 
contributions to the theory of multirate signal processing. 

Another opportunity for further work on multirate systems is in the area of 
lattice structures. Some work on lattice structures has already been researched, but 
extensive development of lattice theory as it applies to multirate linear prediction can 
lead to more convenient representation and implementation. In addition, this disser¬ 
tation focused on optimal solutions using time-domain characteristics. Solutions to 
optimal hltering of multirate systems using frequency-domain techniques would fur¬ 
ther expand the theory and utility of multirate statistical signal processing. Research 
into hlter design of periodically correlated scalar time series using a spectral factor¬ 
ization algorithm for wide-sense stationary vector processes has already been started 
by Spurbeck and Scharf [Ref. 29]. 

Research from a companion thesis [Ref. 50] has extended the concepts of opti¬ 
mal multirate hltering to two-dimensional (2-D) signal processing and high-resolution 
image reconstruction. The reader of this thesis may want to refer to [Ref. 50] for 
some additional topics not treated here. 
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APPENDIX A. POWER SPECTRAL DENSITY 
FOR MULTIRATE PROCESSES 


This appendix contains the derivations of the power spectral density functions 
for the the decimator, expander and hlters presented in Chapter If. The standard 
equation for the power spectral density function of a stationary process 

OO 

S(^) = Rll]e-^“‘ (A.l) 

/ = —OO 

can be applied to the input and output correlations independently. However, any 
difference between sampling rates of the input and output signals makes it impossible 
to develop a power spectral density function for the cross-correlation dependent on 
a lag only. Therefore it is necessary to use a two-dimensional power spectral density 
function. This function is dehned as 

OO OO 

,S2^(Ai,Ao)= (A.2) 

ni=—OO no=—OO 

For the input to a decimator, expander or hlter, the power spectral density is 
dehned as 

OO OO 

Sl^{Xi,Xo)= (A.3) 

ni=—OO no=—OO 

Using the time-lag form (when appropriate), the power spectral density function 
becomes 

OO OO 

Sl^{\u)= Y, (A.4) 

n=—oo l=—oc 

which reduces to (A.l) when the signal is wide sense stationary. 
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A. DECIMATION 


The operation of decimation is depicted in Fig. A.l and defined as 


y[niiy] = x[Lmy\ for m = ..., —1, 0,1, • 


A^x]- 


iL 


{«,] 


Figure A.l. Decimator 

The correlation functions for the decimator are as shown in Table 2.2 of Chap¬ 
ter II. The output power spectral density function for decimation is 


OO OO 


Sf(A,,A„)= ^ ^ 


my=—oo my=—oo 
OO OO 


E E Rx[Lmy, Lm'y 




OO OO 




^-jXim'y^-jXomy 


^-j{XilL)Lm' -j{XolL)Lmy 


my=—oo m'=—oo 


Q'2.U(\ \ \ _ q2D I -^1 -^0 

1^1, '^Oj — Ox i 


In the time-lag representation, the output power spectral density function is dehned 
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as 


oo oo 

Sf{\u;)= E 

my=—OQ ly=—oo 

oo oo 

= Y1 R^[Lmy-Lly\e-^^^ye-^^^y 

ITly - OO ly ^—oo 

OO oo 

= ^ ^ R^[Lmy-,Lly\e-^^^R^^^ye-R‘^R'^^^y 

my = — 00 ly = — 0 O 


s“(A.u,) = si;‘'|p^ 


or in lag representation only 


OO 

S,{u,) = E 

ly - OO 

OO 

= R4Lly]e-^^^y 

Ly - OO 

OO 

= Rx[Lly\e-R^IR^^y 

Ly - oo 

S,(oj) = S. (j) 


The cross-power spectral density function is given by 

OO OO 

Sl?{\x,\y)= J2 Y1 Rxy[m..my]e-R^^^e-Ry^y 

mx=—oo my=—oo 
oo oo 

= ^ ^ R^[m^,Lmy]e-R^'^^e-Ry^y 

mx=—oo my=—oo 

oo oo 

= ^ ^ R,[m,,Lmy]e-R-^-e-R^yR^^^y 

mx=—oo my=—oo 
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In the time-lag representation, the cross-power spectral density function is dehned as 

OO CO 

mx=—oo ly=—oo 
CO OO 

TTlx — — *0O 


OO OO 

^ ^ ^ ^ Rx['^X') Lly 

RRy — *0O —OO 

OO CO 

^ ^ ^ ^ Rx\p^x^ 

my=—oo ly=—oo 


Sl^{\x,UJy) = S\ 


'xy 


2D 

X 



^ , ^y\ 




B. EXPANSION 


The operation of expansion is depicted in Fig. A.2 and dehned as 


y[m] = 


x[m/I] m div I 


otherwise 




t / 




Figure A.2. Expander 


The correlation functions for the expander are as shown in Table 2.3 in Chapter 
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II. The output power spectral density function for expansion is 


OO OO 


^f(Ai,Ao)= Ry[my,m'y]e-^^^^ye-^^^^'y 


my=—oo m'=—oo 


E E 


my=—oo m'=—oo 


OO OO 


Y1 


my=—oo m'=—oo 


^f(Ai,Ao) = S' 2 ^ (/Ai,/Ao) 


In the time-lag representation, the output power spectral density function is dehned 


Sl^{\u)= Y RyV^y'^h 


\e •' “e 


TTXij — CXD Iffj — — OO 




m„=—OO ly=—oo 


E E Rx[^y/1] 


-j(I\)my/I^-j(lLO)ly/I 


(AM 

A lag representation for the expander does not exist, since the output is not wide 
sense stationary. The cross-power spectral density fnnction is given by 


OO OO 


Sly{Xx,Xy)= Y E Rxy[m.,rny]e-^^^^^e-^^y^y 


mx=—oo m-u=—oo 


Y Y Rx[m,,my/I]e-R^^^e-Ry^y 


mx=—oo m-u=—oo 


OO OO 


Y, Y, Rx[mx,my/I]e 


-j>^xmx -j{IXy)my/I 


mx=—oo my=—oo 


Sl?{Xx,\y) = Sr{Xx,^ 
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In the time-lag representation, the cross-power spectral density function is dehned as 

OO CO 

mx=—oo ly=—oo 
CO OO 

TTlx — — *0O 

OO OO 

fR'y — OO —OO 

OO CO 

Rx[m^; {ly + (/ - 

my = — OC ly = —00 

s^^;f{x.,Uy) = sl^ {K -{I- 1K,M) 




C. FILTERS 

Filtering is depicted in Fig. A.3 and dehned as 

OO 

y[m\= h[r\x[m — r\. 


x\m\ 



y[m\ 


Figure A.3. Filter 


The correlation functions for the hlter are as shown in Table 2.4 in Chapter 
II. The output power spectral density function for a hlter is 

OO OO 

^f(Ai,Ao)= X] E Ry[m,rn']e-R^^e-R^^' 

m=—oo m'=—oo 

CO OO 

= {Rx[m, m'] * h[m] * h*[m']) 

m=—oo m'=—OO 

(Ai,Ao) = ^r(Ai,Ao)i/(Ai)i/*(Ao) 
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In the time-lag representation, the output power spectral density function is dehned 


as 

OO OO 

Sf(A,a,)= ^ ^ 

m=—OO /=—OO 
OO OO 

= ^ ^ {Rxlm; 1] * h[m] * h*[—l]) 

m=—oo /=—OO 

or in lag representation only 

OO 

S„{u,) = R,ll]e-^“‘ 

/ = —OO 
OO 

/ = —OO 

Sy{uj) = Sxiu;)H{u)H*{uj) 

The cross-power spectral density function is given by 

OO OO 

S“(Ai,A„)= E 

m=—OO m'=—oo 
OO OO 

m=—OQ m'=—oo 

S“(Ai,A„) = Sy(Ai,A„)/f(A„) 

In the time-lag representation, the cross-power spectral density function is dehned as 

OO OO 

S“(A,a,)= Y1 E 

m=—OO /=—OO 
OO OO 

m=—oo l=—oc 
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D. SUMMARY OF CORRELATION AND POWER SPEC¬ 
TRAL DENSITY 

For ease of reference, a summary of the correlation functions and power spec¬ 
tral densities for decimation, expansion and filtering is provided in Tables A.l through 

A.3. 


Table A.l. Summary of Decimation 




iL 




Ry[my,m'y] = R^[Lmy, Lm'y\ 

sy(A.,A„) = sj®(y,^) 

Fty'^TTly^ /^] Rx\_L7Tly^ vl 

sy(A,^) = sy (A 

^yih] ~ Rx[Lly] 

S,(a;) = S, (1) 

RxylR^xj Rx\R^xj 

S“(A.,A,) = Sy (a„%) 

RxylR^xi ^y\ Rx\R^xi ^^y 

Sxy ^Xx + L 
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Table A.2. Summary of Expansion 


x[m^] 


t / 



II 

R^[my/I, m'y/I] niy, m'y div I 

0 otherwise 

5f(Ai,Ao) = 52^(/Ai,/Ao) 

II 

R^[my/I]ly/I] my, ly dlV / 

] otherwise 

Sl^{\u:) = Sl^{I\,Iu:) 

Rxy\j^X) 

j Rx[mx, my/ 1 ] my div / 

1 0 otherwise 

Sl/;{K,\y) = Sl^{K,IXy) 

Axy [Ula; , /y] = ^ 

Rx[m:^] {ly + (/ - l)m,^)/I] m^ - ly div I 

0 otherwise 

Sl/;{K,U:y) = Sl^ (A, - (/ - l)^y, iLOy) 


Table A.3. Summary of Filtering 
x[m\ -► h[m\ - >■ y[m\ 


Ry[m, m'] = Rx[m, m'] * h[m] * h*[m'] 

(Ai,Ao) = 52^(Ai,Ao)i^(Ai)Tr*(Ao) 

Ry[l] = Ry[l]*h[l]*h*[-l] 

Sy{u) = Sx{uj)H{u;)H*{u) 

Rxy[m, m'] = Rxim, m'] * h*[m'] 

Sl^{\i,\o) = Sl^{X^,\o)H*{Xo) 

Rxy[m, 1] = Rx[m] 1] * h*[-l] 

S“(Ai,A„) = Sy(A,^)ff-H 
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APPENDIX B. KRONECKER PRODUCT AND 
REVERSAL NOTATION 


A. KRONECKER PRODUCT 

This appendix provides the definition of a Kronecker Product and lists some 
useful relations. A thorough discussion of Kronecker products and matrix calculus is 
contained in [Ref. 62] and a shorter concise synopsis is contained in [Ref. 63]. 

Given two matrices A and B, the Kronecker product A 0 B is defined as 

fliiB 012 B ... oijB 

a2iB a22B ... a2jB 

A = 

ttiiB aj2B ... ttijB 

Some basic properties and relations for Kronecker products are provided in Table B.l. 

Table B.l. Some Properties of Kronecker Products 

(1) (A 0 B) (D ® G) = AD (g) BG 

(2) (A + H)(g)(B + R) = A(g)B + A(g)R + H(g)B + H(g)R 

(3) (A ® B)^ = A^ ® B^ 

(4) (A (g) B)“^ = A-i 0 B-i 

(5) tr (A (g) B) = tr (A) ■ tr (B) 

(6) (A (g) B) (g C = A (g) (B (g) C) 
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B. REVERSAL NOTATION 


Given a vector v of length P whose elements are 

r iT 


Ol O2 


CLp 


the reversal of the vector, denoted by v, is dehned by 

r iT 


V = 


(2 p 


ap_i 


cti 


Given a P x Q matrix A with elements aij, 




Ol,l 

Ol,2 ■ ■ ' 

■ Oi,Q 

A 

= 

02,1 

02,2 ■ ■ ' 

■ 02,Q 



Op,l 

Op,2 ■ ■ ' 

■ Op,Q 

the reversal of A is given by 






Op,Q 

Op,Q-l 



r Op-1,Q Op-1,Q-i ■ ■ ■ (ip- 1,1 

J\. — 

Oi,Q oi^Q-i ■ ■ ■ ai,i 

A few important properties of vector and matrix reversal are listed in Table B.2. 
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Table B.2. Properties of Reversal (from [Ref. 18]) 



Quantity 

Reversal 

Matrix product 

AB 

AB 

Matrix inverse 

A-i 

(A)-i 

Matrix conjugate 

A* 

(A)* 

Matrix transpose 

A^ 

(A)^ 
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APPENDIX C. DERIVATION OF 
PARAMETERS FOR INNOVATIONS FORM OF 
THE OPTIMAL FILTER 


In Chapter III recursive optimal filters were used to present the optimal hlter 
in its innovations representation. This appendix contains the derivations related to 
these optimal hlters. 


A. REPRESENTING THE OPTIMAL FILTER IN A RE¬ 
CURSIVE FORM 


The Wiener-Hopf equations pertaining to an optimal multirate hlter with M 
input channels are given by (3.7). The solution for the hlter coefficients can be written 


as 


where 


and 









■p (k) 

Kii 

■p (k) 


p (k) 
^MM 


'^d(M+l) 




In general, for the signal (C.2) and (C.3) can be dehned as 


(C.l) 


(C.2) 


(C.3) 





(C.4) 


R 


(k)*T 


R 


(fc) 
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and 


^di — 


;(fc) 

-dl 


;(fc) 

■ d{i—l) 


In order to derive a recursive form of the optimal filter, the inverse 
is first written in a recursive form. By partitioning into 





g(t, 

f,(k)*T 

rW 


where 


= 


R 


(fc) 

ij 


■p(fc) 


the recursive form of (see [Ref. 69]) can be expressed as 


= 


/~i(k) ^ 
^i-1 

0 

+ 

0 

0 


'^j-1 

0 

+ 

0 

0 



-G. 

I 


(fc) 


E 


(fc)- 


_ ^{k)*T j 


(fe)l^(fe) ^ ^(k) *T _^ 


Cj. "^E- G 


I t 

^-l 


E 


(fc)- 


where 




and 


eS‘> = RSf> - Bif'cfV bS*' 

^ii • 


(C.5) 

matrix of 

(C.6) 

(C.7) 

(C.8) 

(C.9) 
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Applying the result of (C.7) specifcally to the signal by substituting (C.7) 
into (C.l) leads to 


K 


(fc) 


,(fc) 


0 0 






■plfc) 


Separating the two components of then produces 


{k)*T 


E 


(fc)- 

M 







^dM 



Ak) 



'-dM 


Ak) 

'■M 


p(fc) ^ Q 

0 0 






'-'dM 

Ak) 

+ 


_ dM _ 



(fc)-|-^(fc) ^ f~i{k)*T f~^{k)T^{k) 




M 


'M 


E 


(fc)- 

M 





"dM 


Ak) 


_ dM _ 


Extracting the terms associated with , the filter is found to be equiv¬ 


alent to 


u(^) _ _ p(^) ^p(^)*^u(^) _1_ p(^) 

~ '^M ’^dM ^ ‘■dM 


_ p(^) ^ ( z.ik) p(^)*^C.(^) 

— [‘■dM '^M ^dM 


(C.IO) 


and the rest of the matrices reduce to 


.(fc) 


Ak) 

^M-1 


_ p(^) _i_ 

~ '^M-1 ’^dM “T 


p(^)p(^) ^(~^ik)*T 


p(^)p(^) ^ 


Ak) 

■^dM 

Ak) 

-dM 


or 



h 


(fc) 

M-l 


pl«; Ul«J 1 p 
'^M-1 ^dM + 


pV^V ^p 
M 


■'(fc) _i_ pi(^)p(fc) ^ f^ik) *T-r{k) _ p'-'^Jp 

'^dM 


(fc)p(fc) ^~(fc) 


M 


M 


-dM- 


Further manipulation leads to 



h 


(k) 

M-l 


p(fc) 

'^M-1 ^dM 


p (^)p(^) 



-G 


(k)*Tk(k) \ 
M ^dM j ’ 


(C.ll) 
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and substituting (C.8) and (C.IO) into (C.ll) produces 


(fc) 


hi 


,(fc) 


_ 

~ '^M-1 '^dM ^M-1^■ 


(C.12) 


Equation (C.12) can also be expressed as 


h(fc) 


_ p(fc) ^ 




p-1 

'^M-1 


I'd(M-l) 

and substituting (C.7) into (C.13) results in 


b(^) 

b(^) 


h(fc) 


(C.13) 


h(fc) 


'^M-2 

0 0 


+ 


(fc) 

'^M-1 


■p(fc) ^p((fc)*r 


'M 


xr(^) ^ 


(fc) *r 
'M-l 


-1 


e(^) 
^M-l 







^d(M-l) 



^d{M-l) 


r^{k) ^ Q 
'^M-2 

0 0 


+ 


'^M-l^M-1 '^M-1 


(^) ^ 
''^M-l^M-1 


^M-1 '^M-1 


{k)*T 


Extracting the terms associated with leads to 

_P’(^) ^ 


E 


-1 


(fc) 
M-l 


B 


(k) 


.(fc) 


(M-l)M^'-M 


h(fc) _ 


■^M-l^d(M-l) 


p(fc)-ip(fc)*ri3(fc) 
^M-l'^M-l ^(M- 


u(fc) 


E(fc)’ *j?' 
^ M -1 ■‘^1 


(fc) u(fc) 

{M-1)M'-^M 


This equation can be expressed equivalently as 

u(fc) _ ( f. 

“M-1 ~ ^M-l 1 ^ 


:(fc) 

d{M-l) 


p(fc)*rr(fc) 


_ Tp(fc) ^ / R >- 

^M-1 \^{M-1)M 


,{k) 


(fc) *T-ry{k) 

'^M-1 ^(M-1)M I ‘■‘■M 


(k) 


(C.14) 
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The terms for through 


are then written as 



h 


(fc) 

M-2 


c 


M-2 


(k) 

d{M-l) 


+ 


'^M-1 


(^) ^ 
'^M-l^M-1 


'(k) 

^d{M-l) 

™(fc) 


^M-2^(M-1)M^‘-M 


(k) 




B 


(fc) 


hfc) 


{M-l)M‘-'~M 


This equation can be rearranged, producing 


hS^) 

uik) 

^'■M-2 


c 


M-2 ^ 


(fc) 

d{M-l) 


'^M-l^M-1 


(f{k) 

y-d{M-l) 


-G 


ik)*Tr{k) 

M-1 ^d(M-l) 


+ G 


(fc) 

M-l 


(fc)- 

'M- 


R 


(k) 

(M-1)M 


(k) ^r-pj (fc) 
'^M-1 ^(M-l)M 



- c 


(fc) 


-1 


M-2 


h 


(fc) 

M • 


(C.15) 


Now substituting (C.8) and (C.14) into (C.15) yields 


h) 


(fc) 


h(fc) 

^‘■M-2 


“ '-'M-2 ^d(M-l) 


'-'M-2 ^{M-l){M-iy'-M-l 


f^{k) ^■d(^) 
^M-2^{M-1)M'-'-M ■ 


(C.16) 


By examining the pattern of how (C.l) is separated into (C.IO) and (C.12) and how 
(C.12) is separated into (C.14) and (C.16), the following recursions can be observed 




E 


(fc)- 


E 


(fc)- 

M 


U{k) 

y-di 

- Gl-’-^bW) 

M 


Y. 

Ef " (R|f - 

j=i+i 


{'■dM 

p.(fc)*rr(fc) 

'^M ’^dM 


qfc)*Tg(fc)\ 


(1 < i < M - 1) 


(C.17) 

(C.18) 
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where the matrices and are given by 


Gf ^ = i 


and 


■p (k) 

Kii 

p {k) 

p {k) 


p(fc)*r 

■'^12 


p(fc) 

^2{*-l) 



fl(k)*T 
^2{j-l) ■ ■ 

1 

7 


Ef ) = 

1 rSI 

i = 

1 


-1 


i = 1 


Br i<i<M 


Rjf - Gf ^ 1 < i < M 


The hlter in (C.17) and (C.18) can be fnrther simplihed by dehning 


(fc) _ -pCfc)"' fftik) _ Q{k)*T-Q{k) 


W = E. 


fik) 

‘■di 




G 


(fc)*Tr(fc) 
i '^di 


(C.19) 

(C.20) 


and snbstitnting these terms into (C.17) and (C.18) 

M 

hf) = hf^ (1 < i < M - 1) (C.21) 


j=i+l 

h£' = h'^f (C.22) 

The hlter coefficients of (C.21) for an M-signal system can be expressed in 
terms of the cross terms and orthogonal hlters only, which is usefnl for 
development of the recnrsive error variance. Beginning with the general eqnation for 
the hlters for signals 1 to M — 1, (C.21), bnt nsing the index ii instead of j for the 
snmmation 

M 

h^^ = hY+ Y (-Hj^)hSf (1<^<M-1) (C.23) 

*i=*+i 


the hlter for signal i can be expanded by snbstitnting (C.21) for the hlter h. 


(fc) 

*1 


Ak) 


A{k) 


M 


+ E 

*i=*+i 


H 


*1 


M 


+ E -h.7 


(k) \ ^(fc) 

"*i*2 ; *2 


12=*!+! 


(1 < i < M - 1). 

(C.24) 
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Separating the expanded term yields 

M MM 

hl‘)=h:<‘'+ 5; (-Hlf) hy + 5; ^ (-Hj') (-H® ) h<f (l<i<M-l). 

il=*+l ii=i+l i2=*l+l 

(C.25) 

Substituting (C.21) into the third term leads to 

M 

hf’ = h;'‘'+ (-<')h;f 

*i=*+i 

MM f M \ 

+ ^ 5: (-hS‘')(-H®) h;f'+ 5: (-HS:>)h<‘' (l<«<Af-l) 

*1=*+1 *2=*1 + 1 V *3=*2 + l / 

(C.26) 

Separating the last term on the right-hand side of the equation yields 

M MM 

hf> = h;'->+ 5: (-<>)h;«+ 5: (_H«) (-H.«) hif> 

il=i+l ii=i+l *2=*l-|-l 

M M M 

+ E li E (-<’) (-Hit’.) (-Hlti) hit (1 < i < M - 1) (C^ 27 ) 

11=1+1 i2=*l+l *3=*2 + l 

If this recursion is continued, the hnal result becomes 

M MM 

hf> = h;'->+ 5: (-<>)h;«+ 5: 5: + ... 

*l=*+l il=*+l 12=*!+! 

MM M M 

y ... y y V■ ■ 

Z^ Z^ y “v V / V *M-2*M-1 ) y *M-1*M ) *M 

il=i+l *2=*1+1 iM-l=*M-2 + l *M=*M-1+1 

(l<i<M-l). (C.28) 

The result of (C.28) can be used to also develop the recursive form of the error 
variance. 


B. EXPRESSION FOR THE ERROR VARIANCE 

Beginning with the basic dehnition of the error variance 

M 

<^I = -Rj(0)-EbEh!‘’*. (C.29) 
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the summation can be separated into 


M-l 

”1 = ^ hif, (C.30) 

i=l 

Substituting (C.28) into (C.30) for the filter term representing the hlters for signals 
1 through M — 1 leads to 


M-l 


al = i?rf(0) - f. 


(k)T 

di 


X 


h. 


i=l 

M MM 

'(fc) , f(k) , 

iii 

\ 

il=i-\-l 22 = 21+1 


21 = 2+1 


z (-hs?) + 5: E (-hS‘ 

il=i+l i2=ii+l 

MM M M 

, 5 : E ■■■, z , E 

ii=i+l *2=*i+l *m-i=*m-2+1 *m=*m-i+1 






-H® I X 


H 


(k) 

* 1*2 


H 


(k) 


*M-2*M-1 


H 


(fe) 




h. 


'(fc) 


_ (fc)* 

'-dM '■^M ■ 


(C.31) 
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The multiple summation term can be separated into two parts leading to 


M-l 


al = i?rf(0) - f. 


(k)T 

di 


X 


i=l 



M-l 



*1=1+1 


M-l M-l 


E E (-h£')( 

* 1=*+1 * 2 =* 1+1 




+ ■ ■■ 


M-l M-l M-l M-l 

+ E E ■+ E . E 

*1=*+1 *2=*1+1 *M-1=*M-2 + 1 *M=*M-1+1 



+ E E-, E 


-H 


(k) 

'* 1*2 


H 


(k) 


*M-2*M-1 


1 h 


'(fc) 


~{k)T.{k)* 

MM ^^M ■ 


(C.32) 


The first two terms on the right-hand side is the error variance for an (M — 1)- 
signal system. Thus dehning crl = al and representing the error variance for the 
(M — l)-signal system as al leads to 


a 


2 

k,M 


= a 


2 

k,M-l 



{k}* 

M 


M-l 



*=1 




M-l M 

-E E 47(-<>•) 

i=\ ii=i-\-l 

M-l MM M 

-EE.E-. E 47(-h.T)x 

i=l ii=i-\-l ^ 2=^1 + ! 2 + 1 



H 


^M-2^M-1 


] h 

"*M-iM 1 


'{k)* 


(C.33) 
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C. SUMMARY 


For clarity and ease of reference, the results for the innovations form of the 
optimal hlter are summarized here. Given an M-signal multirate system, the solution 
to the Wiener-Hopf equations is given by 


hf)' 


rG) 

rvii 

■p (k) 

-1 

1 - 

<—1 

1 



p (^) 

■ ^MM 


z.ik) 

‘■dM 


The following terms are dehned for use in the hlter equations 


= 


i 

bG) _ 

1- 

-■ 

1 

_1 

^di — 

-1 


and 


and 


Gf ^ = ( 


p (^) 
Kii 

p {k) 
■•^12 

p(fc)*r 

■*^12 





Ef’ = 


■p (^) 

rtii 


i = 1 


R 


r: 


(fc) 

ip-i) 

(fc) 

2{i-l) 


R 


(fc) 


bP l<i< M 


i = 1 


rI^^ - of ^ l<i<M 


The equations for the hlter coefficient in the innovations representation of the hlter 
are then summarized in Table C.l. 
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Table C.l. Innovations Form of the Optimal Filter 
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APPENDIX D. MULTICHANNEL LEVINSON 

RECURSION 


This appendix briefly discusses the procedure known as the multichannel 
Levinson algorithm, (also called the Levinson-Wiggins-Robinson (LWR) algorithm). 
For more detailed discussions of this topic see [Ref. 64, 70, 71]. The multichannel 
Levinson algorithm is a method of solving the Normal equations for multichannel sys¬ 
tems in a recursive fashion. This algorithm is an extension of the Levinson recursion 
for single-channel systems. 

A. THE MULTICHANNEL SYSTEM 

In order to properly describe the multichannel Levinson algorithm, it is nec¬ 
essary to first consider a multichannel random process and its correlation function. 
Let the random process x[n] have the form 

xi[n] 

X2[n 

xc[n]\ 

where C is the number of channels. The correlation function associated with x[n] is 
given by 

Rx[^] = T{x[n]x*^[n — /]}. 
where Rx[/] is a C x C matrix. 

The multichannel Levinson algorithm solves for the Alter coefficients and pre¬ 
diction error covariance matrices for linear prediction of the multichannel random 


x[n = 
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process. In so doing, it finds the parameters for both the forward and backward 
prediction of the multichannel random process. 

B. GENERAL FORM OF THE ALGORITHM 

To initialize the recursion, the following variables are dehned, 

R-O = Rx[l] Ag = Bg = IcxC Sg = Eg = Rx [O] , 

where Icxc is an identity matrix of size C x C. In addition, a cumulative correlation 
matrix is dehned 

Rx[p] 

Rx[p -1] 

Rx[l] 

where p is between 1 and the desired order. Throughout this algorithm, the variables 
with the superscript b refer to the backward prediction parameters and the other 
variables refer to the forward prediction parameters. 

Step 1. Calculate the rehection coefficients, Tp and T^ 

r 1 ~ rr i ~ 

Ap = [Rf ... r^Ja;_, = [[r, r2 ... RpJb;_, 

Tp = Sj:\Ap 

and 

■pb _ p-1 A T 

Step 2. Calculate the order prediction coefficients, Ap and Bp 




148 





Bp_i 0 

K= _ - _ K 

0 

Step 3. Calculate the order Error Variance coefficients, Sp and 



These steps are repeated until the desired order has been reached. At that 
point, the prediction coefficients, error variances and reflection coefficients for all 
orders up through the desired order have been calculated. 
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APPENDIX E. A GENERALIZED MATRIX 
TRIANGULAR FACTORIZATION 
ALGORITHM 


In order to solve (4.41) for the transform matrix Aq when developing a so¬ 
lution to the multirate linear prediction problem in Chapter IV, a generalized UL 
factorization was developed. The generalized UL factorization presented in this ap¬ 
pendix is modeled after the standard LU factorization algorithm that can be found in 
many linear algebra books, such as [Ref. 65]; however, when factored, the left matrix 
is upper triangular and the right matrix is lower triangular. The algorithm is called 
“generalized” since it allows for a block UL factorization using different size blocks. 
It thus provides for Gaussian elimination of multiple rows at the same time. This is 
important to the multirate linear prediction solution, since eliminating multiple rows 
at the same time is necessary when dealing with joint observations of multiple signals. 

A. LU TRIANGULAR FACTORIZATION REVIEW 

This section briefly reviews the algorithm for computing the basic LU trian¬ 
gular factorization of a matrix. For a more detailed discussion of this methodology, 
refer to any elementary linear algebra textbook, e.g, [Ref. 65]. 

Given an m x n matrix A, 

flip fli,2 ■ ■ ■ fll,n 

fl2,l 0-2,2 ■ ■ ■ 02,n 

Om,l fl *,2 ■ ■ ■ Om,n 

it is desired to decompose the matrix into upper and lower triangular matrices, such 
that 

A = LU, 
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where L is the lower triangular matrix with I’s on the main diagonal and U is the 
upper triangular matrix. The procedure for this method is as follows: 

Step 1: For i = 1,2,.. . ,m — 1, eliminate row i from rows i + 1 to m: 

a. Calculate row multipliers for rows i + 1 to m 

Mk,i = ak,ilai^i , k = i + l,i + 2,... ,m 

b. Subtract row i from rows i + 1 to m 

ak,i = ak,i - Mk,i * ai^i = 0, k = i + l,i + 2,... ,m 
ak,i = ak,i - Mk,i * a*,/ , k,l = i + l,i + 2,... ,m 

After eliminating i rows, the matrix has the form; 

flip fll,2 fll,3 ■ ■ ■ fll,i fll,i+l ■ ■ ■ fll,n 

0 02,2 fl2,3 ■ ■ ■ 0-2,i fl2,i+l ' ' ' fl2,n 

0 0 03,3 ■ ■ ■ fl3,i fl3,i+l ' ' ' fl3,n 

A, = ■ ■ ■ ■ ■ 

b b b • • • 0^,2 O^,2-1-1 * * * 02,21 

b b 0 ■ ■ ■ 02+1,2 fli+l,i + l ■ ■ ■ fll+l,Tl 

b b 0 ■ ■ ■ flm,i flm,i+l ' ' ' flm,n 

Repeat step 1 until matrix A is upper triangular. This is the matrix U. 


fll,l fll,2 fll,3 ■ ■ ■ fll,i fll, 2 +l ■ ■ ■ fll,n 

b 02,2 fl2,3 ■ ■ ■ 02,i 02,2+1 ' ' ' fl2,n 

b b 03,3 ■ ■ ■ 03,2 03,2+1 ■ ■ ■ 03,21 

U = A,2i = 

b b 0 • • • Oi,! fll,l + l ' ' ' fll,21 

b b 0 ■ ■ ■ 0 02+1,2+ 1 ■ ■ ■ fli+1,21 

b b 0 ■ ■ ■ 0 0 ■ ■ ■ 0271,21 
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Step 2: The lower triangular matrix can be formed as 


L 


1 0 

M2,i 1 


0 

0 

1 


0 

0 

0 , 


Mm,2 Mm,3 ' ' ' 1 

where Mij are dehned above. 


B. A GENERALIZED UL TRIANGULAR FAGTORIZA- 
TION ALGORITHM 

The basic LU triangular factorization algorithm can be modified and extended 
to a generalized UL triangular factorization algorithm. This new procedure requires 
the dehning and use of interim variables to transform an original matrix A. 

Given an m x n matrix A, 

0 - 1,2 ■ ■ ■ ai,n 

02,2 ■ ■ ■ 02,n 

it is desired to decompose the matrix into generalized upper and lower triangular 
matrices, such that 

A = UL, 

where U is a block upper triangular matrix and L is a block lower triangular matrix. 
The blocks are not all of the same size, but the sizes of the desired blocks are known 
(or given). 

The procedure for this method is as follows (starting with the last column), 


A = 


Ol,l 

02,1 

Om,l 
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Step 1: Eliminate columns i — r + 1 through i from columns 1 through i — r, where r 
is the number of columns to be eliminated at the same time and i is the index of the 
last of these columns. 

a. Dehne the following vectors and matrices 


r+1,2—r+l 


A 










r+l,fc 




, A;= 1,2, 


^i,k 




r+1 



1 , 2 , 




— r 


b. Calculate column multipliers for columns 1 to i — r 

M{i-r+l)i,k = ^ — 1, 2, . . . , f — r 

c. Subtract columns i — r + 1 through i from columnss 1 through i — r 


^{i-r+l)i,k ^{i-r+l)i,k ■^{i-r+l)i,{i-r+l)iM(^i-r+l)i,k 0 , k 1, 2, . . . , f T 

(^l,k ®Z,fc ^l,(i—r+l)i-^{i—r+l)i,k ) k^ I 1, 2, . . . , f T' 


At the iteration, the matrix has the form: 
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^1,2—r 

■ ■ T 

+ 


0'l,n-2 

^l,n—1 

<^1,22 

®3,1 

Cti —— f- 



Oji—r,n—2 

Oj—r+l,n—2 

^2—r,n—1 

^2—r+l,n—1 

^i—r,n 

^2—r + 1,22 












®i,n—2 

^2,22—1 


0 

0 

0 

0 ■■ 

Clm—2,m—2 

2 , 22—1 

^222 — 2,22 

0 

0 

0 

0 ■■ 

0 

^222—1,22—1 

^222—1,22 

0 

0 

0 

0 ■■ 

0 

0 

^m.n 


Repeat step 1 until matrix A is transformed into the upper triangular matrix U. 




^1,2 

fll,3 ■ ■ 

Ol,i 

Ol,i+l 

*^1,22 


0 

0-2,2 

02,3 ■ ■ 

02,i 

02,1+1 

*^2,22 


0 

0 

03,3 ■ ■ 

03,i 

03,i+l 

*^3,22 

U = A^ = 

0 

0 

0 ■■ 

Oi,i 

Oj,i+l 

^i,n 


0 

0 

0 ■■ 

■ 0 

Oj+l,i+l 

*^2+1,22 


0 

0 

0 ■■ 

■ 0 

0 

^m,n 


Step 2: After the 2'^'^ row is eliminated from the 1®*^ row, the upper and lower trian¬ 
gular matrices can be formed as 


L 


1 0 0 

1V[2^1 1 0 

M3^2 1 


0 

0 

0 


Mm,2 Mm,3 ' ' ' 1 
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C. EXAMPLE 


Given a matrix A, 


A = 


flip 

: ai^2 

fll,3 • 

: 01^4 : 

■ fll,5 

fll,6 

fl2,l 

: fl2,2 

fl2,3 : 

: 02,4 : 

■ fl2,5 

fl2,6 

fl3,l 

: fl3,2 

fl3,3 : 

: 03,4 : 

: fl3,5 

fl3,6 

fl4,l 

fl4,2 

fl4,3 • 

! CI4 4 ! 

: fl4,5 

fl4,6 

fl5,l 

fl5,2 

fl5,3 

fl5,4 : 

: fl5,5 

fl5,6 

fl6,l 

fl6,2 

fl6,3 

fl6,4 : 

: fl6,5 

^6,6 


(E.l) 


where ajj is the matrix element for the row and column. It is desired to decom¬ 
pose the matrix using the generalized UL algorithm. For this particular example, it 
also is desired to eliminate columns 5 and 6 together and followed by the elimination 
of columns 2 and 3 together, as shown by the partitioning in (E.l). The number of 


12 12 


, where 


columns to be eliminated can be represented in vector form as 
the numbers refer to the block sizes (number of columns to be eliminated together). 
Step 1: Eliminate columns 5 and 6 

a. Dehne the following vectors and matrices 
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^56,56 


0-5,5 ^5,6 

06,5 06,6 



05,1 


05,2 


^5,3 


^5,4 

^56,1 — 

^ 56,2 — 

^56,3 — 

^56,4 — 


06,1 


06,2 


^6,3 


06,4 


^1,56 — 


Ol,5 Oifi 


^2,56 — 


02,5 02,6 


^3,56 — 


^3,5 ^3,6 


^4,56 — 


^4,5 Oa6 


b. The multipliers for columns 1 through 4 are 


Mse,! — M5g^2 — A5g^5ga5g^2 

M56,3 = Agg^ggagg^s Mgg^4 = Agg^ggagg^4 

c. Eliminating columns 5 and 6 from columns 1 through 4 produces (primes 
have been added to the modihed values, for this example only, to highlight the changes 
at each step) 


! 

■56,1 ~ 

^56,1 — 

A56,56M56,1 = 


= ^1,1 

— ai,ggMgg,i 

®3,1 

= 03,1 

~ a3,ggMgg,i 

/ 

■56,2 = 

^56,2 — 

A56,56M56,2 = 

! 

Oi 2 

= ^1,2 

— ai,ggMgg,2 

/ 

®3,2 

= ^3,2 

~ a3,ggMgg,2 

/ 

■56,3 ^ 

^56,3 — 

A56,56M56,3 = 

®1,3 

= fll,3 

— ai,ggMgg,3 

! 

®3,3 

= 03,3 

— a3,ggMgg,3 

/ 

■56,4 ^ 

^56,4 — 

A56,56M56,4 = 

! 

Oi 4 

= Ol,4 

— ai,ggMgg,4 

/ 

®3,4 

= 03,4 

~ a3,ggMgg,4 


®2,1 ~ 02,1 — a2,56Mgg^l 

04 1 = 04,1 — a4,ggMgg,i 

®2,2 ~ ®2,2 — a2,56Mgg,2 

04 1 = 04,2 — a4,ggMgg,2 

02,3 = 02,3 — a2,56Mgg,3 

®4,3 ~ ®4,3 ~ ^4,56M56,3 

<^2,4 = 02,4: — a2,56Mgg,4 

<^4,4 = ^4,4 — a4,ggMgg,4 
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d. Matrix A with columns 5 and 6 eliminated looks like 



. / 

■ ®1,2 

®1,3 

. ! 

■ ®1,4 

: ®1,5 

fll,6 

/ 

®2,1 

. / 

• ®2,2 

/ 

®2,3 

. ! 

• ®2,4 

• ®2,5 

0 - 2,6 

®3,1 

. / 

■ ®3,2 

! 

®3,3 

. / 

■ ®3,4 

: 03,5 

03,6 

/ 

®4,1 

/ 

®4,2 

! 

®4,3 

. / 

. ^ 

: ^4,5 

04,6 

0 

0 

0 

0 

: 05,5 

05,6 

0 

0 

0 

0 

: ^6,5 

06,6 


Step 2: Eliminate row 4 

a. The multipliers for columns 1 through 3 are 

l = Cl^ \l ^4^4 -Tf4^2 ~ ^ 4 , 2 /®4,4 -^4,3 ~ ^ 4 , 3 /®4,4 


b. Eliminating column 4 from columns 1 through 3 produces 


// 

®4,1 = 

/ 

®4,1 ~ 

®4^4-Tf4^1 — 

// 

®2,1 

= ®2,1 

02,4'^^4,1 

// 

Oa,2 — 

/ 

®4,2 ~ 

O^^^Mi,2 — 

// 

02,2 

/ 

— 02,2 

(^2, 4 -^ 4 ,2 

// 

®4,3 ~ 

/ 

®4,3 ~ 

Oy^y^MA,3 = 

// 

®2,3 

/ 

~ ®2,3 

— 02^Ma,3 


// / ' 1\ /T 

^3 1 — ^3 1 ^3 4 '^^ 4,1 

Oj^ ^ CL^ 4-^^4 1 

// / ' 1\ yT 

®3,2 ~ ®3,2 ~ ®3,4^4,2 

// / ' 1\ yT 

Qj^ 2 ~ 2 4 -^ 4 ,2 

®3,3 ~ ®3,3 ~ ®3,4^44^3 

(Xi 3 — *^1,3 *^1^4-^4,3 
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c. Matrix A with columns 4 through 6 eliminated looks like 


tt • tt tt • t 


4,1 ' 

: 2 

®1,3 ■ 

: 0^4 : 

• ®1,5 

Ol,6 

// 

■ // 

// 

. / 



'2,1 ' 

: 02 2 

®2,3 ' 

: 02 4 : 

• ®2,5 

02,6 

// 

■ // 

// 

. / 



b,i ' 

' ®3,2 

®3,3 ' 

: O34 : 

: 03,5 

03,6 

0 

0 

0 : 

/ 

CI4 4 

: 04,5 

04,6 

0 

0 

0 

0 ! 

: 05,5 

05,6 

0 

0 

0 

0 

: 06,5 

06,6 


Step 3: Eliminate columns 2 and 3 

a. Dehne the following vectors and matrices 


^23,23 — 

// 

02,2 

// 

// 

02,3 

// 

^23,1 — 

// 

02,1 

// 

^1,23 — 

// // 
Ol,2 O43 


03,2 

03,3 


03,1 




b. The multiplier for column 1 is 

M23,1 = ^23^23^23,1 

c. Eliminating columns 2 and 3 from column 1 produces 


^23,1 


— ^23,1 — A23,23M23,1 — 0 


1,1 


= a 


1,1 


ai,23M 


23,1 
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d. Matrix A with columns 2 and 3 eliminated looks like 


/// 

. // 

: a-^^2 

// 

®1,3 

. / 

■ ®1,4 

: ®1,5 

Ol,6 

0 

• // 

• ®2,2 

// 

®2,3 

. / 

• ®2,4 

• ®2,5 

02,6 

0 

• // 

■ ®3,2 

// 

®3,3 

. / 

■ ®3,4 

: 03,5 

03,6 

0 

0 

0 

• / 

• ^4 4 

: ^4,5 

04,6 

0 

0 

0 

0 

: 05,5 

05,6 

0 

0 

0 

0 

: 0-6,5 

*^6,6 


Step 4: Build the U and L matrices 

a. The matrix U is the transformed matrix A 


/// 

Oil : 

. // 

• ®1,2 

// 

®1,3 

. ! 

■ ®1,4 

: Ol,5 

Ol,6 

0 

• // 

: 02 2 

// 

®2,3 

. ! 

• ®2,4 

■ ®2,5 

02,6 

0 

• // 

• ®3,2 

// 

®3,3 

. ! 

■ ®3,4 

: 03,5 

03,6 

0 

0 

0 

• ! 

. Ct^ ^ 

■ ^4,5 

04,6 

0 

0 

0 

0 

: 05,5 

05,6 

0 

0 

0 

0 

■ ^6,5 

^6,6 
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b. The matrix L is made up of the multipliers, such that 

r 1 ; 0 0 : 0 : 0 0 

1 0 : 0 ; 0 0 

0 1 : 0 ; 0 0 

Mii^2 Tf4^3 : 1 : 0 0 

; 1 0 

M56^2 M56^4 

; 0 1 


M23,1 

L = ... 
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